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ABSTRACT 


Low frequency combustion instability, known as chugging, is 
consistently experienced during shutdown in the fuel and oxidizer 
preburners of the Space Shuttle Main Engines. Such problems always 
occur during the helium purge of the residual oxidizer from the 
preburner manifolds during the shutdown sequence. Possible causes 
and triggering mechanisms are analyzed and details in modeling the 
fuel preburner chug are presented in this thesis. A linearized 
chugging model, based on the foundation of previous models, capable of 
predicting the chug occurrence is discussed and the predicted results 
are presented and compared to experimental work performed by NASA. 
Sensitivity parameters such as chamber pressure, fuel and oxidizer, 
temperatures and the effective bulk modulus of the liquid oxidizer are 
considered in analyzing the fuel preburner chug. The computer 
program CHUGTEST is utilized to generate the stability boundary for 
each sensitivity study and the region for stable operation is 
identified. 
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CHAPTER I 


INTRODUCTION 


The Space Transportation System (STS), better known as the Space 
Shuttle, is propelled into orbit by two solid rocket boosters and 
three main engines (Fig. 1.1). The Space Shuttle Main Engines (SSME) 
have successfully completed twenty-four flights and several hundred 
test stand firings. The propellants utilized by these engines are 
cryogenic liquid hydrogen and liquid oxygen. The engines are stable 
for steady- state operation as well as programmed load changes from 
minimum power level (MPL) to full power level (FPL). There have, 
however, been some problems associated with the SSME during shutdown 
both in flight and on the test stand. In the shutdown sequence, the 
engine oxidizer system is purged by heliiim prior to fuel cutoff. 

During the helium purge, the engines experience a low amplitude low 
frequency pressure pulsation (chugging) due to combustion instability 
in the fuel and oxidizer preburners. Since the thrust is essentially 
reduced to zero there have been no significant effects on the engine's 
performance. These pulsations have been linked to undesirable 
bearing loads and damage to the augmented spark ignitor (ASI) line. 
This problem may require frequent replacement of engine components and 
can potentially make the STS less cost effective. Therefore, to reduce 
cost, NASA is sponsoring this research in an effort to reduce the 
maintenance costs of future SSME designs. 
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Figure 1.1 Space Transporation System Schematic. 

Source: Sutton G. P. , Rocket Propulsion Elements , John Wiley and 

Sons, New York, 1986. 







Space Shuttle Main Engines Description 

Fig. 1.2 is a schematic of the propellant flow in the SSME show- 
ing the major equipment. Two preburners are used on each engine to 
preheat the hydrogen and supply power for the fuel and oxidizer pumps. 

Cryogenic liquid hydrogen enters the engine from the external 
tank via a low pressure prunp which supplies enough head to prevent 
cavitation of the three stage high pressure fuel pump. Following the 
high pressure fuel pump the hydrogen leaves at approximately 7000 
psia, at full power level, and is used to cool the nozzle, throat and 
main combustion chamber before entering the preburner at a temperature 
of 1027 K. The majority of the fuel flows through the preburners and 
is partially burned before entering the main combustion chamber". 

The oxidizer (liquid oxygen) follows a similar path but is not 
used for cooling purposes and enters the main combustion chamber 
directly. The pressure at the outlet of the single stage pump of the 
oxidizer preburner is approximately 4400 psia and is fed directly into 
the main combustion chamber operating at 3277 psia. A portion of the 
oxidizer is supplied to the prebumer p\mp which supplies oxygen at 
6850 psia to the fuel and oxidizer preburners. 

The main engines supply a nominal thrust of 470,000 Ibf at rated 
power level. Engine power is controlled by throttling the oxidizer 
flow to the two preburners via the fuel and oxidizer preburner 
oxidizer valves (FPOV and OPOV). The power available to the 
turbopvimps and hence reactant flow rate is controlled by the supply of 
oxygen to preburners 
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.2 Schematic of propellant flow in the Space Shuttle Main 
Engine. 

Sutton G. P. , Rocket Propulsion Elements . John Wiley and 
Sons, New York, 1986. 


Source: 




Under steady-state conditions the preburners are operating at 
very fuel rich conditions with an approximate equivalence ratio of 
six. Combustion is extinguished in the preburner and in the main 
combustion chamber by closing the OPOV in the oxidizer line. The 
closing of the OPOV causes the suspension of the oxidizer flow to the 
combustion chamber thus extinguishing combustion prior to fuel cutoff. 

The engines are throttled back prior to shutdown in flight 
operations due to maximum acceleration limitations of the shuttle. 

The ground test firings on the SSME are frequently shutdown from one 
hundred percent rated power level. Fig. 1.3 shows the shutdown 
sequence for a typical engine operation near one hundred percent of 
rated power level. The oxidizer preburner oxidizer valve (OPOV) is 
closed first. Shortly after the closure of the OPOVj the fuel 
preburner oxidizer valve (FPOV) closes. The preburners are isolated 
from the oxidizer system once the two preburner oxidizer valves are 
closed. The only oxidizer available to the preburners is the 
residual trapped in the line and manifold volvime between the valve and 
combustion chamber. This residual oxygen is cleared into the 
preburners by a heliiim purge. 

Fig. l.A shows the purge and the augmented spark ignitor piping 
for the fuel preburner. A similar arrangement exists for the 
oxidizer preburner. The check valve in the helium purge supply shown 
in Fig. 1.4 remains closed until the pressure downstream of the valve 
drops to about 750 psia. The residual oxidizer is cleared from the 
ASI line, the oxidizer valve and line and the oxidizer manifold into 
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the preburner combustion chamber where combustion with the fuel takes 
place. Although the fuel flow rate varies due to pump speed varia- 
tion, substantial fuel flow is maintained until after the purge is 
completed. 

The pressure pulsation experienced during the helium purge of the 
oxidizer in the fuel preburner is called chugging based on the rela- 
tively low frequency (75--200 Hz) and its apparent cause which is a 
coupling between propellant feed rate and combustion generated 
pressure in the preburner. Chugging usually begins about 2.3 to 2.5 
seconds after the cutoff command is given on ground test. The condi- 
tions for flight test data are not available; however, the situation 
should be similar to ground testing except that in flight cut-off 
usually occurs at lower power levels. It may appear from Fig. 1.3 
that the chug start correspond to the main oxidizer valve (MOV) 
closing, however, the preburners are completely isolated from the 
oxidizer system by the time MOV closes. 

Problem Description and Method of Solution 

The SSME consistently experience a low frequency pressure pulsa- 
tion, called chugging, in both the fuel and oxidizer preburner combus- 
tion chambers during the helium purge following cut-off. The area of 
interest is shown in Fig. 1.5. The amplitude, frequency and duration 
of the instability appear to be dependent on shutdown conditions 
particularly the helium compressibility and fuel temperatures. The 
chug occurs, only on shutdown of the SSME, during the helium purge 
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Figure 1. 
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of the fuel preburner oxidizer manifold. 

Since the chug occurs on shutdown, loss of performance on the 
SSME is of no major concern. The chug peak pressure never exceeds 
1000 psia in the combustion chamber designed for pressures in excess 
of 6000 psia. Because of this fact, chugging has received very 
little attention. However, failure of the augmented spark ignitor 
line and turbopump bearings has been linked to the chug. 

The primary objective of this research is to improve on analysis 
techniques for liquid propellant rocket instability in the fuel 
preburner of the SSME. The triggering mechanism and the possibili- 
ties of elimination were also studied. System variables such as 
chamber pressure, fuel flow rate, oxidizer flow rate, and fuel and 
oxidizer temperatures were varied and their effect on the fuel 
prebumer chug is presented in the following chapters. 

Combustion instabilities, particularly in rocket engines, have 
been the subject of several previous investigations. The analysis 
presented by Crocco and Cheng [14] and Harrje and Reardon [7] are more 
extensive in the analysis of combustion instabilities compared to 
earlier models. Their models concentrated on design criteria that 
should be avoided during engine construction. The early developments 
of combustion instability analysis were linearized models based upon 
the work done by Crocco and Cheng. These models have shown that 
chugging is critically sensitive to combustion time delay and low 
injector pressure drops. The compressibility .of the feed system was 
also considered as a principal contributor to the chug. A literature 
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search by the author showed that most chugging problems are analyzed 
using the one-dimensional Itunped parameter model. The approach taken 
by J. R. Szuch [12] in analyzing main combustion chamber chugging 
instabilities appeared to be the most effective of the techniques 
available. The computer program written by Szuch was used as a 
foundation for this work. The method was adapted to the SSME 
propellant flow Arrangement including the unchoked exit flow in the 
SSME fuel preburner. 

Following the analysis performed by Szuch [12] the characteristic 
equation, derived from the non-linear differential equation modeling 
the propellant flow, describing the stability of a bipropellant rocket 
system is reduced to a quadratic equation. The quadratic formula, 
which is a closed-form solution, is used to solve the characteristic 
equation. A high accuracy solution is obtained as no iterative 
process is required. The two roots obtained are either real and 
distinct, equal or complex conjugates. The solution of the characte- 
ristic equation, utilizing critical system values, will result in the 
generation of a critical stability boundary. It is therefore 
possible to predict if the operating point of the fuel preburner is 
in the stable or unstable region. Only the complex roots with posi- 
tive real parts are sought. Complex roots with positive real parts 
indicate that the system is unstable; the imaginary part obtained 
provides the frequency of interest at the critical operating point. 

A chug analysis program was written for the fuel preburner of the 
SSME and implemented on the VAX 11/780 computer at the University of 
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Tennessee. The program provides insight into the critical parameters 
affecting chugging and suggests ways to avoid operating conditions 
conducive to chugging . The code was verified by comparison to 
several experimental test cases which are discussed below. Program 
results in the form of stability diagrams will be analyzed and 
compared to data obtained from static firings performed under the 
direction of NASA engineers at Rocketdyne and National Space 
Technology Laboratory. 
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CHAPTER II 


LITERATURE REVIEW 


A literature review of liquid propellant rocket instability was 
conducted by the author to gain some insight into the problems and 
work associated with combustion instabilities. Although the majority 
of this research is devoted to low frequency instability analysis, a 
brief discussion on the different types of rocket instability is 
presented. 

Instability refers to an uncontrolled pressure oscillation in the 
rocket combustion chamber. In extreme cases, characterized by an 
increase in amplitude with time, the combustion chamber may rupture or 
be severely damaged. Rocket instabilities are characterized by the 
frequency of the pressure oscillation. There are basically three 
types of rocket instabilities; low, intermediate and high frequency 
referred to as chug, buzz and scream respectively. These 
instabilities result from the coupling between the combustion of the 
propellant and the fluid dynamics of the system. Low frequency 
instability involves a coupling between the combustion chamber 
pressure and the propellant flow rate. Intermediate frequency 
instability involves standing pressure waves in the feed line system. 
High frequency instability involves the occurrence of pressure waves 
in the combustion chamber. 
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One such type of instability, low frequency instability, which 
occurs during shutdown of the Space Shuttle Main Engines (SSME) is the 
principal focus to be discussed in this thesis. Chugging is 
characterized by the constructive coupling between the combustion 
energy release and the propellant feed process. The chugging 
frequency is generally less than 200 Hz, although there is no clear 
distinction in the cutoff frequency range between low and intermediate 
combustion instabilities. 

Fig. 2.1 shows the schematic of a mono -propellant rocket system. 



Figure 2.1 Schematic of mono-propellant rocket system. 

Assuming that the injector flow velocity (v^) decreases by a small 
amount, this will cause a decrease in the chamber pressure (F^) which 


14 



will increase the flow velocity. Shortly after the input flow 
increases, will also increase thereby decreasing v^. This process 
continues and the pressure oscillation is known as chugging. The 
amplitude could increase or decrease with time. If the amplitude 
increases with time, the system is unstable and continued chugging 
will cause severe damage and undesirable irregularity in the thrust of 
the rocket motor. 

Most of the research on chugging was done in the early 1950 's 
with the birth of the space program. The chugging phenomenon was 
observed by Summerfield [1] during a series of tests on a 1000 Ibf 
thrust rocket motor. For mathematical simplicity, Summerfield 
analyzed the mono -propellant rocket system neglecting the inertia of 
the liquid propellant. Utilizing the concept of Crocco's [2] time 
lag theory, the governing equation was derived from the conservation 
of momentum and energy. The resulting non-linear differential 
equation was linearized by perturbation methods and expressed as: 

u” 4- A u’ + Bu + Ct’ (f — c) = 0 (2.1) 


where v is the flow acceleration of the propellant, r is the combus- 
tion time lag and A, B and C are parameters evaluated at the operating 
conditions. The general solution of eq. (2.1) can be represented as 
a linear sum of the particular solutions as: 


y {A. +• jGi ) t 

n = 0 


( 2 . 2 ) 
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The roots of the characteristic equation ( X , ) are obtained 

by substituting eq. (2.2) into eq. (2.1), grouping real and imaginary 
terms, and transforming into: 


2 2 — X 

X +ar + 6 — qr + ce cos q; = 0 


(Real) (2.3) 


(2.x + a)qi — ce ^ sin q;= 0 


(Imaginary) (2. A) 


where a, b, c, x and ^ , are system parameters expressed in terms of 
the combustion time lag. 

Plotting these equations on the x and V plane, shown in Fig. 
2.2, it is possible to locate the intersections of equations (2.3 and 
2.4). The behavior of the the rocket system can be determined for 
any given values of A, B, C and 'P . The solutions involving 
positive values of x lead to instability; while negative values of x 
represent stability. 



Figure 2.2 Approximate representation of equations (2.3 and 2.4). 

Source: Summerf ield, M., "A Theory of Unstable Combustion in 

Liquid Propellant," ARS 21, 1951, pp. 108-114. 
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It is the objective of the designer to select appropriate values 
of A, B, and C such that these values are damped out. Therefore, it 
is possible to manipulate values of A, B, and C such_that solutions 
with positive values of x are impossible. This leads to a sufficient 
condition for stability defined by Summerfield [1] as: 


m 

P A„ 

c 2 


+ 


C*L* 
R T 

g c 


2AP 


> T 


(2.5) 


* ■ * 

where C and L are the characteristic velocity and length 
respectively, 1^ and m are the length of the feed line and propellant 
flow and A? is the pressure drop across the injector. 

From eq. (2.5) the following actions may be inferred for over- 
coming instability: 

1. increase the pressure difference between the supply tank and 

combustion chamber either by reducing the area of injector orifices or 

inserting resistance elements in feed lines 

* 

2. increase L , which is the ratio of combustion chamber volume 
to throat area 

3. increase the length of the feed lines to the combustion 
chamber 

4. reduce the cross-sectional area of the propellant feed line 

5. reduce the combustion time lag ( r ) of the propellant either 
by changing to a more reactive propellant or adding a catalyst to the 
propellant 
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Summerf ield [1] noted that changing the propellant from nitric 
acid--gasoline to nitric acid--aniline removed the instability. The 
aniline, being self-igniting reduces the combustion time delay, con- 
firms that different propellants should be examined. 

The analysis performed by Gunder and Friant [6] suggests ways to 
eliminate the chugging instability of a bipropellant rocket motor. 

The governing differential equation was derived from the equation of 
motion of the propellants. Four different methods of solving the 
characteristic equation were proposed. The Nyquist stability 
method, which involves a conformal transformation, was chosen. 

Gunder and Friant [6] showed the instability could be determined, 
by plotting the transformation of the governing differential equation 
in the f(z) -plane. If the plot encircles the origin the system is 
unstable. In the plots shown below. Fig. 2.3 shows an unstable 
system, while the system shown in Fig. 2.4 is stable. The system 
represented in Fig. 2.4 is similar to that in Fig. 2.3 but was stabi- 
lized by changing certain parameters in the system, such as, 
decreasing the area of feed lines or increasing the pressure drop 
across the tank and chamber for the oxidizer and fuel. 

Several authors including Summerfield and Gunder and Friant have 
approached the chugging instability of liquid propellant rocket with a 
simplified mathematical description of the combustion process to 
obtain a closed-form solution or a simplified differential equation 
for analog and digital computers. Webber [4] presents a method of 
calculating low frequency unsteady combustor behavior based upon the 
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Figure 2.3 Nyquist diagram for unstable rocket system. 



Figure 2.4 Nyquist diagram for stable rocket system. 

Source: Gunder, D. F., and Friant, D. R., "Stability of Flow in 

A Rocket Motor," Journal of Applied Mechanics, 17, 1950, pp. 327-333. 
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simultaneous integration of the differential equations describing: 

1. the propellant flow in the feed system and the injectors and 

2. the evaporation and combustion rates of propellant droplets 
The model included a fuel and oxidizer supply system, an injector, a 
combustion chamber and a nozzle. The calculations are based upon the 
Euler marching technique for the calculated rate of change of . the 
system variables and the integrating time interval. 

In Webber's model, the time-varying injection rate of each 
propellant component is expressed as a function of chamber pressure. 
The lumped parameter approximation approach is used to obtain: 


dt 


AP - l-Rpq" 
pS//A 


( 2 . 6 ) 


for the liquid propellant flow incorporating only resistance and 
inertia terms, where 1 and A are the length and cross-sectional area 
of the feed lines, R is a constant to approximate line losses and q 
is the volxometric flow rate.- The summation is performed only if the 
feed line has varying cross-sectional area or multiple feed lines are 
used. The initial axial velocity of the droplet produced by the 
injected stream is; 


V =K2u. 

I j 


(2.7) 


where the constant K2 includes the velocity losses due to atomization 
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( 2 . 8 ) 


process and is the velocity of the jet spray expressed as: 

m 

j pA 

The atomization process is approximated by computing the time-varying 
mean diameter as a function of physical propellant properties and 
system geometry. Webber calculates the mean diameter of the droplet 
(D^) empirically from: 


D 

m 


7.99 [^yJptD. 
39.37 + 195 D.} 


(2.9) 


where is the injector orifice diameter in meters, p, 5 and P 
are viscosity, surface tension and density respectively. The 
respective units are defined in the list of symbols. 

The evaporation rate of a stable propellant is dependent upon the 
heat transfer mechanism to the surface of the droplet and the 
molecular diffusion away from it. The heat transfer mechanism from 
the surface to the interior of the droplet is significant, however, 
the assiomption that heat transfer to the interior of the droplet is 
infinitely fast, due to small droplet size, simplifies the analysis. 
From this assumption the rate of evaporation of the propellant droplet 
is: 

lil = A'u nk D ^ (2.10) 


where Nu is ,the Nusselt number and dH is the sum of sensible heat 



and latent heat required to raise the droplet temperature from the 
injected temperature to the boiling point, and k is the thermal 
conductivity of the gas. The rate of efflux of gas from the chamber 
is expressed in terms of time-varying pressure, temperature and ther- 
mal properties in the combustion chamber. The exhaust gases are 
assxjmed to be ideal, well mixed and in thermo chemical equlibrium. 

Two interesting results reported by Webber [4] are shown in 
Figures 2.5 and 2.6 below. Both plots were calculated for a 100 msec 
sequence of high-amplitude, 66 Hz, chugging instability. The 
integrating time interval used was 0.25 msec. Webber discovered that 
the numerical calculations appear to be negligibly influenced by the 
integrating time interval for the nominal values used. When a double 
time interval of 0.5 msec was used the amplitude changed by only 2 
percent and the frequency changed by less than 3 percent. The pro- 
gram was tested in comparison with an experimental engine producing 
low frequency instability. The experimental and calculated 
frequencies were found to differ by not more than 7 percent with a 5 
percent agreement in amplitude. 

Few researchers have analyzed chugging instability using 
different vaporization rates for the propellant. Wenzel and Szuch 
[8] postulated a model with different delay times (vaporization times) 
applied to the respective propellants. The difference between this 
model and previous models lies in the feed system. In the previous 
analysis, the dead times for all the process were Iximped together and 
applied equally to both propellants. While it is true that mixing 
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Figure 2.6 Variation in fuel flow rate during chugging. 

Source; Webber, W. T., "Calculation of Low Frequency Unsteady 
Behavior of Liquid Rocket from Droplet Combustion Parameters," ARS 26 
1956, pp. 26-39. 


and reaction times are common, vaporization times associated with the 
individual propellants should be treated separately. The feed system 
in this model is assumed to be completely decoupled from the combus- 
tion chamber. Figures 2.7 and 2.8 show the comparison between the 
two models. 

The governing equation is derived from the conservation of mass 
and the Laplace transformation is used to obtain: 


0 = 1 + 
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( 2 . 11 ) 


Assuming that the burned gas behaves ideally with constant gas 
residence time ( 6g ) and incorporating further mathematical reduc- 
tions, the characteristic equation is reduced to: 




6 S-r- 1 
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al 

iF 



£ + 1 ? C* 

g 



AP. Ie -I '■« C' 
n g 


( 2 . 12 ) 


for implementation on a digital computer solving for engine stability. 
Here is the gas mixing time, and are the vaporization time of 
the respective propellants, is the mixture ratio and al and bl are 
constants evaluated from steady (mean) operating conditions. 

This model indicates that stability may be achieved by 
manipulating the vaporization time, shown in Fig. 2.9, from 0.5 to 2.0 
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Figure 2.7 Block diagram for single vaporization rate model. 
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vaporization 
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Figure 2.8 Block diagram for proposed model using different 
vaporization rates. 

Source: Wenzel, L. M., and Szuch, J. R., "Analysis of Chugging 

in Liquid Bipropellant Rocket Engines using Propellants with different 
Vaporization Rates," NASA TN — D--3080, 1965. 
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Figure 2.9 Effect of vaporization time on stability boundary. 

Source; Wenzel, L. M., and Szuch, J. R., "Analysis of Chugging 
in Liquid Bipropellant Rocket Engines using Propellants with different 
Vaporization Rates," NASA TN — D — 3080, 1965. 

msec at a constant mixing time of 1.0 msec. If the operating condi- 
tion is as indicated, then increasing the oxidant vaporization time 
from some lower value up to 1.5 msec improved the stability. This 
behavior is contrary to that predicted by the single delay time model, 
since an increase in vaporization time is destabilizing. Further 

increase in vaporization time, however, shifts the curve upward resul- 
ting in unstable behavior. The stability was achieved due to the 
relatively short mixing time. As the mixing time was increased to 
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4.0 msec, increasing the oxidant vaporization time was destabilizing. 
There are limited experimental data to substantiate the behavior 
predicted by this model, however, this model should not be abandoned 
as more experimental data will becomes available. 

Barrere and Moutet, [3] did experimental analyses on low frequency 
combustion instability with two objectives in mind: 

1. experimental determination of the influence of various 
parameters on an engine with low frequency instability, and 

2. comparison of theoretical and experimental results 

They performed a series of experiments involving both a rectangular 
and a circular combustion chamber. Their analysis showed that the 
nature of the propellant, characteristic length of the combustion 
chamber and the mean pressure of the flame body are important para- 
meters in influencing the chugging instability. High injection 
pressure, mixture ratio or length of combustion chamber did not 
influence the instability significantly. The experimental results 
obtained were consistent with Crocco's [2] theory, as long as the 
chamber pressure takes on the same value as Crocco's theory. The 
results were inconsistent with Crocco's theory as soon as the chamber 
pressure was altered. 

Another attempt at eliminating combustion instability lies in the 
modification of the fuel supply line proposed by Li [5]. A suitable 
servo-mechanism is designed to control the fuel flow which is contro- 
lled by the pressure oscillation in the combustion chamber. In his 
analysis, Li treats the dynamic performance of the fuel supply system 
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as a first order system assximing an incompressible fuel supply. The 
purpose of using the feedback system is to associate familiar techni- 
ques to stabilize the feedback system such as that shown in Figures 
2.10 and 2.11. 




Combustion 

Chamber 

P 

rU6X 

Supply 


Stabilizer 
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Figure 2.10 Schematic diagram for Rocket System 



Figure 2.11 Block diagram for Rocket System 

Two methods are presented to improve the stability of the first 
order system. The first is to decrease the time constant so that at 
the unstable frequency the phase lag was reduced. The other method 
is to increase the time constant so that at the iinstable frequency the 
amplitude ratio will be reduced. There is an undesirable increasing 
phase lag associated with the second approach due to the increase in 
time constant. 

Two types of stabilizers were proposed: 

1. flow rate actuated-type stabilizer, and 
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2 . 


flow acceleration actuated stabilizer 


The construction of the flow rate actuated stabilizer, shown in 
Fig. 2.12, uses a spring-loaded piston valve and a sensing nozzle to 
regulate the propellant flow. A pressure drop across the nozzle will 
produce a force, on the piston valve, that will open or close the 
valve port. The flow acceleration actuated stabilizer, shown in Fig. 
2.13, uses a spring-loaded piston valve to control two sets of ports. 
The piston has an annular cavity which allows the entire propellant 
flow to pass through. The inlet and outlet ports are arranged such 
that the propellant flows perpendicular to the axial motion of the 
piston. The movement of the piston is thus due to the change in 
momentum of the propellant in the annular cavity. 

The solution to the chugging instabilities with this approach is 
not as simple as it seems. There are problems associated with the 
operation of both stabilizers at the chugging frequency. Furthermore, 
this theory has not been verified experimentally. 

Some of the ideas presented by Harrje and Reardon [7] to 
eliminate low frequency instabilities include: 

1. increasing the pressure drop in the injector 

2. increasing the fluid inertance (i.e. longer 1/d ratio in the 
injector or feed system) 

3. decreasing chamber volume 

A. changes in the time delay in the combustion of the fuel, 
such as vigorously mixing the recirculated hot gases with the incoming 
propellant to preheat the fuel, and 
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Source: Li, Y. T., "Stabilization of Low Frequency Oscillations 

of Liquid Propellant Rocket with Fuel Line Stabilizer," ARS 26, Jan. 
1956, pp. 26-39. 
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5. 


increasing the damping process to dissipate the oscillatory 
energy in the combustion chamber or decreasing the coupling between 
the driving forces. 

The idea of decreasing the chamber volume, which seems to contra- 
dict Summerf ield' s work, is because Harrje and Reardon define the 
combustion time delay of the fuel as: 


c = 


I 


(2.13) 


where 


t 
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yR t 

^ 8 




LU. j2(Y- 1) 


(2.14) 


Therefore decreasing the chamber volume is similar to decreasing L 
and increasing c and in essence agrees with Summerf ield' s theory. 

The model presented by J. R. Szuch [12] showed promise for 
analyzing the fuel preburner chug with minor modifications. In this 
model the governing equation is obtained from the conservation of mass 
applied to the combustion chamber with the ideal gas assiimption and 
sonic nozzle exit. The model was tested with varying fuel injector 
annulus area, analogous to injector pressure drop changes, and the 
results obtained were compared to experimental and analog data. The 
results obtained from the computer model compared very well with the 
analog data and are presented in Fig. 2.14 below. 

In the literature review the author has discovered that most of 
the work and research done on rocket instability analysis involved 
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Figure 2.14 Stability boundary generated by analog and digital 
computers . 

Source; Szuch, J. R. , "Digital Computer Program for Analysis of 
Chugging Instabilities," NASA TN--D--7026, Dec. 1970. 


linearized steady-state models. This approach is applicable to the 
SSME fuel preburner chug, however, it does not give any information on 
the magnitude of the chug amplitude which are transient and non- 
* linear. The model presented by Szuch [12] was selected and 
applicable for the SSME fuel preburner analysis. The mathematical 
formulation and details of the analysis are presented in the following 
chapter . 
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CHAPTER III 


FUEL PREBURNER CHUG MODELING 


Analysis 

The model presented by Szuch [12] for the analysis of combustion 
instabilities has been modified and is now capable of analyzing chug- 
ging instabilities in the fuel preburner of the Space Shuttle Main 
Engines (SSME). This chapter is thus devoted to the mathematical 
formulation used for modeling the fuel preburner chug. A satisfac- 
tory working model, capable of predicting fuel preburner pressure 
excursion, is presented in this chapter. Although this model 
analyzes the fuel preburner chug, it could also be applied to the 
oxidizer preburner with minor modifications. 

Although detailed analysis of the actual combustion process 
taking place in the fuel preburner is beyond the scope of the present 
work a few elementary assvimptions allow analytical treatment. The 
burnt gases in the fuel preburner can be considered as an ideal gas, 
and the flow is determined by the conservation laws of mass, momentxim 
and energy. In the range of low frequencies pertaining to chugging, 
the propagation of the pressure waves is assumed to be instantaneous, 
thus, the following fundamental assumptions for chugging can be 
stated: 

1. the gas pressure is uniform throughout the combustion chamber 
and oscillates about a mean or steady-state value 
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2. the temperature of the gases is constant and uniform 
regardless of the pressure oscillations, but depends on stoichiometry 

3. the time lag (or delay) is uniform, that is, it has the same 
value for all the propellant elements, and 

A. the combustion of the propellant is infinitely fast once the 
droplet has evaporated ^nd mixed, therefore chemical kinetics may be 
neglected. 

These four assvimptions greatly simplify the analytical treatment 
of the combustion process; the introduction of assumptions (l) and 
(2) replace the momentiim and energy equations which therefore need not 
be considered. The introduction of assumptions (3) and (A) allow 
the finite droplet evaporation time to collapse to a single 
characteristic time. With the combustion proceeding infinitely fast, 
a quasi -steady state process is assumed. 

These assumptions reduce the formulation for the dynamics of the 
combustion process in the fuel preburner to be essentially governed by 
the balance of mass, that is, the rate of burned gas produced must be 
equal to the sum of the rate of ejection of gas out of the fuel 
preburner through the high pressure fuel turbine and the rate of mass 
accumulated in the preburner. 

Derivation of the Characteristic Equation 

The combustion process of the fuel preburner during shutdown in 
the SSME is shown schematically in Fig. 3.1. Following shutdown, the 
oxidizer line is closed by a ball valve and the residual oxidizer is 
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Fig. 3.1 Schematic of Fuel Preburner Combustion 

purged by the helium flow. The fuel and helium flow rates are deter- 
mined by conditions upstream and downstream of the injector elements. 
The oxidizer vaporizes, mixes and reacts with a small fraction of the 
fuel to produce a hot gas consisting primarily of and H 2 O. The 
time interval that exists between liquid propellant injection and 
conversion to vaporized propellant is referred to as the vaporization 
time delay. The characteristic vaporization time delay ( ), 

used in this model, is the time required for the conversion of 50 
percent of the liquid propellant to the vaporized state. The gas 
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phase mixing time delay ( a ) is the time interval required for con- 
version of vaporized propellant to burned gases. These time delays 

( o and o ), which are fxinctions of time, are assumed to be combined 
mu 

into a single characteristic time. The following equations, there- 
fore, relate the rate of change of burned products to the injected 
propellant flow rates. 


m AO = m .{t — a — o ) 

00 Qt VO m 


(3.1) 


m„ iO — m..At — a . — a ) 

fb fi vf m 


(3.2) 


where m is the mass flow rate, and subscripts ob, fb, oi, and fi 
represent oxidizer burned, fuel burned, oxidizer injected and fuel 
injected respectively. 

Using a one-dimensional lumped parameter approach, Harrje and 
Reardon [7] state that the basic governing equation is obtained from 
the conservation of mass. The conservation of mass written for the 
fuel preburner assviming that all the reactants are burned yields: 


£ 

dt 



= m-(0 + m ,(0 — m (0 
fo 00 e 


(3.3) 


where p is the gas density, is the volume of the fuel preburner 
combustion chamber; subscript e represents the exit flow. . The helium 
purge is included in the m^j^ term. Eq. (3.3) is based on the assump- 
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tion that the volume occupied by the liquid is very small compared to 
the total chamber volume, and is valid only for low frequencies [7], 
The exit mass flow is given by [9]: 


m = CP 

e c 



1.4286 



1.7143 

IpJ 



(3.4) 


where C is an empirically determined coefficient which is a function 
of turbine speed; P and T are the pressure and temperature 
respectively. Subscripts c and hg represent the combustion chamber 
and hot gas manifold respectively. 

Since the gas density is a ftmction of pressure and temperature, 
eq. (3.3) is non-linear and must be linearized to obtain the desired 
solution. Linearization of eq. (3.3) is performed by assuming a 
small perturbation in the system variables about a steady-state 
operating point. Neglecting the products of perturbations 
yields the following equation: 
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(3.5) 


where the curl superscripts are the perturbation quantities. 


0 = 
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(3.6) 
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is the gas residence time, and 


m = m,4- m m 
t fb 00 i 


(3.7) 
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is the total mean mass flow of the exhaust gases. The process 
involved in acquiring eq. (3.5) from eq. (3.3) is found in Appendix A. 
Assuming that the gas behaves as an ideal gas, the first term on the 
left hand side of eq. (3.5) can be shown to be: 
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(3.8) 


Since the temperature perturbation varies, that is, the 
perturbation at the combustion front and the exit of the preburner is 
not the same, the mean temperature perturbation should be integrated 
over the entire preburner combustion chamber. However, with 
assrumption (2) stated earlier the average chamber temperature remains 
relatively constant. Thus, the second term on the right hand side of 
eq. (3.8) goes to zero. 

Utilizing the ideal gas law in eq. (3. A), the second term on the 
left hand side of eq. (3.5) may be expressed as: 


rh W 


= Cl 
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t 




(3.9) 


where Cl is determined from input parameters, and is a function of the 
specific heat ratio of the combustion products and pressure ratio 
found in eq. (3. A). 
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Equations (3.5), (3.8), and (3.9) are combined to yield the 


following equation: 
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(3.10) 


Assuming that is the adiabatic flame temperature then clearly. 




(3.11) 


for any given set of reactants where 4> is the equivalence ratio. 
Therefore, the last term in eq. (3.10) becomes; 
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since the temperature is a weak f\inction of chamber pressure, the last 
term in eq. (3.12) is neglected. The equivalence ratio perturbation 
may be expressed in terms of the oxidizer and fuel mass flux perturba- 
tions as; 
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(3.13) 
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which is analogous to eq. (5.2. 1-7) in Harrje and Reardon [7]. 

Equations (3.1), (3.2), (3.10), (3.12) and (3.13) are combined 
yielding the following equation which relates the chamber pressure to 
the injected propellant flow rate. 
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where t = (o + o )» is the combustion time delay, and 
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o and f refer to oxidizer and fuel respectively. The Laplace 
transformation of eq. (3.14) is; 


9 s + Cl 


P (s) = 

c 


’ XA ,(s)e + 

C Ol 


P Y m As) e 

C ft 


-I J5 
! 


(3.17) 


The perturbation in the injected propellant flow rate is related 
to the perturbation in the chamber pressure at the injector face by: 
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(3.18) 
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(3.19) 


where Z^(s) and Z^(s) are defined as the linearized output impedance 
of the oxidizer and fuel feed system respectively, and is the 
chamber pressure at the injector face. The determination of Z^(s) 
and Zj(s) will be presented in the next section. 

The total and static chamber pressure at the injector face are 
related by their definition as: 


P 


ct 


= Pr 


I + 


- 1 9 

M~ 

2 



(3.20) 


where M is the Mach number and Y is the specific heat ratio. Sub- 
scripts. t and s represent the total and static, pressures respectively. 

The static pressure in the combustion chamber and at the injector 
face are written in terms of the Mach number as: 


p I + Y ‘W" 
Pc l+Y^/^ 


(3.21) 
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Assuming that M . is small, (i.e. a large chamber), equations (3.20) 
cx 

and (3.21) can be combined to yield the following pressure ratio (K) . 
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(3.22) 


Assuming that there is no loss in the total chamber pressure, equa- 
tions (3.17), (3.18), (3.19) and (3.22) are combined to yield the 
following result: 
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The above equation relates the steady- state average chamber 
pressure to the perturbed chamber pressure in the system. The 
existence of the chug is inferred by an unbounded pressure 
perturbation, that is, P^,(s) approaches infinity. This condition is 
obtained by setting the denominator of eq. (3.23) to zero resulting 
in; 


-V 

" + 


-1 = 




Z . i.s) 


9 s -h Cl 


(3.24) 
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Eq. (3.24) is, therefore, the characteristic equation describing 
the stability of the bipropellant rocket system. Eq. 3.24 applies to 
the fuel preburner system when the proper choice of K, P^, X, Y, Z, 
and Cl are made. 

Eq. (3.24) is solved only for complex roots having positive real 
parts (exponential growth). Appearance of these roots indicates that 
chugging is present and defines the stability boundary. The appear- 
ance of a negative root signifies that the system is stable. The 
imaginary root defines the frequency at which the chug occurs. The 
stability boundary (Fig. 3.2) generated for P^ = 4.4815E6 Pa (650 
psia) with T^ = T^ = 40 K (72 R) shows the operating point in the 
unstable region at a frequency of 121 Hz. A minimum frequency of 75 
Hz was selected and the stability program solved the characteristic 
equation for positive real complex roots. If no non-negative roots 
are found, the procedure is repeated with another selected frequency. 
The positive real part of the complex root obtained is related to the 
respective pressure drops ( curve, 

therefore, signifies the boundary at which the positive real part 
exists. The region to the right of the boundary, outside the enve- 
lope, corresponds to stable operation, while the region within the 
envelope is unstable. The shape of the stability boundary is such 
that an increase in oxidizer or fuel pressure drop stabilizes the 
preburner. This type of modeling provides insights into the proba- 
bility of instabilities, but provides no inf oinmation on the limits of 
the chug amplitude. 
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Fig. 3.2 Stability Boundary for Fuel Preburner. 


Solution of Characteristic Equation 


The solution of the characteristic equation is performed by 
separating the output impedances Z^(s) and Z^(s) into their real and 
imaginary parts. Letting s = j co , the characteristic equation is 
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separated, resulting in the following real and imaginary equations 
(3.25 and 3.26) respectively. 
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where 


Z (/co) = R +jl 

(3.27) 

Z^Ajw) = R^ +jl^ 

(3.28) 


For the frequency range of interest, equations (3.25) and (3.26) are 
solved for two critical parameters and R^. The real parts of R^ and 
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are selected because they can be related to the injector flow 
resistances which are independent of the frequency. 

It can be shown that equations (3.25) and (3.26) reduce to a 
quadratic equation relating the critical value of the real part of the 
oxidizer impedance for a specified critical frequency, to the 

imaginary feed system impedances and The procedure, 

outlined in Appendix A, for the combination and reduction of equations 
(3.25) and (3.26), will yield the following quadratic equation: 
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As shown in Appendix A, the solution of eq. (3.29) will be used to 
compute the critical real part of the fuel impedance R^' via: 
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(3.31) 
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The solution of equations (3.29) and (3.31) will be related 
physically to the injector element and feed system being studied. 

The imaginary variables and 1^ must be computed, at each specified 
critical frequency, prior to the solution of equations (3.29) and 
(3.31). The modeling of the feed system, injector element, and 
other system parameters as well as the computation of the operating 
point variables will be presented in the next section. 


System Parameters Computation 

A Fortran 77 computer program based on a program for the 
analysis of chugging instabilities by Szuch [12] was written for the 
fuel prebumer of the SSME. The quadratic formula is applicable for 
the solution of the characteristic equation. The two roots procured 
will either be real and distinct, equal, or complex conjugates. 
Although an iterative solution is not required, the computer program 
is written to calculate the system variables at each specified 
frequency of interest. Thus, a critical frequency is selected and 
the characteristic equation is solved for R^’ and R^’, where the 
superscripts (') represent critical values. The solutions of R^' and 
R^' are then related to the critical oxidizer and fuel pressure drop 
ratios. The main function of the program is, therefore, to generate 
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the stability boundary of the chug during the shutdown transient of 
the SSME for different parametric test cases. Regardless of the type 
of study, the solutions of equations (3.29) and (3.31) require prior 
calculations of the system variables at each operating point. 

The variables X and Y are computed directly using equations 
(3.15) and (3.16) respectively. The steady-state value of exit flow 
rate m^ is the total flow rate of the exhaust gases leaving the 
fuel preburner. Assviming that complete combustion is taking place in 
the fuel preburner, the following reaction equation is applicable to 
the combustion process; 

6 H 2 + 1/2 0^ =======> 5H^ + H^O (3.32) 

The NASA SP-273 code by Gordon and McBride [10] is used, with an 
equivalence ratio ( ) of six, to calculate the combustion chamber 

temperature (T^) . The chamber temperature calculated is 
approximately 1065 K which is too low a temperature to dissociate 
and H^O; hence the assumption of complete combustion is valid. The 
NASA code is also used to obtain a plot of T^ verses. ({> ,. with v^irying 
equivalence ratio from one to eight at the shutdown condition. A 
least squares curve fit is utilized to acquire a polynominal relation- 
ship between T and (p . The slope (^^ / 3 <J)) at the specified 
equivalence ratio can thus be calculated. 

The chamber pressure ratio (K) is computed using eq. (3.22). The 
specific heat at constant pressure (C^) for each component of the 
products is procured from 4 > , and the approximated empirical 
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polynomial expressed as [11]: 


Cph q = 34.19 - 43.87 0°-^^ + 19.78 0°-^ - 0.88 0 (3.33) 

Cp„^ = 13.51 - 167.96 0~‘’-'^® + 278.44 0"^ - 134.01 0'^‘^ (3.34) 

Utilizing the mole fraction weighted average the specific heat ratio 
and the exit Mach nvimber (M^) can be solved from equations (3.35) and 
(3.36) respectively. 


Y = 



(3.35) 


m 

(3.36) 

where is the universal gas constant and the density is obtained 
from the weighted average. 

The gas residence time is computed from its definition (eq. 3.6). 
The vaporization time delay ( ) for the oxidizer is defined as the 

time required to vaporize 50 percent of the mass of the injected 
droplet. The average droplet velocity over this length is 
approximately equal to the injection velocity of the droplet as 
reported by Priem and Heidmann [13]. Based on this information, the 
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vaporization time delay can be calculated from: 



V. 

ip 


(3.37) 


where I^q is the length required to vaporize 50 percent of the 
injected mass droplet, and v^^ is the injection velocity of the 
propellant droplet. The vaporization time delay was computed for the 
liquid oxidizer only. Since the fuel is already completely vapor iz.ed 
before injection into the combustion chamber, the vaporization is set 
to zero. 

The length required to vaporize 50 percent of the injected mass 
is computed utilizing: 


‘50 


T 1 

r 

% |°-"5 

10.35 

[. H ] 

7.62xl0~® 

130.5 1 

100. ) 

3.26x10® J 


0.0699 — (3.38) 

10.66 

2.07 xlO®^ 


where r is the mean droplet radius in meters, M and H are the molecu- 
lar weight and heat of vaporization with units defined in kg/kg mole 
and J/kg respectively [12]. 

The drop radius of the liquid oxidizer is required in computing 
l^Q. Experimental evidence has shown that mean oxidizer drop radius 
for a concentric element, with gaseous fuel injected through an annu- 
lus (shown in Fig 3.3) is proportional to the square-root of the injec- 
tion momentum ratio with a proportionality constant of 0.118 [12]. 
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m V 
o o 


0.5 


(3.39) 


r=K d 

r a 




where is the oxidizer injector diameter. The drop radius is also 

affected by the properties of the oxidizer at the boiling point 
temperature (90 K). This effect is represented by: 


r<x 



(3.40) 


where p is the viscosity and 5 is the surface tension at 1 atm, 90 
K. To calculate the drop radius, must be modified to incorporate 
this effect. Therefore, with the modified value of K^, the resulting 
expression; 


r= O.llSd 

o 


^^oPo2 

0.25 

r m V 

O 0 

^2^02^0 


\mv\ 


(3.41) 


can be used for any liquid oxidizer. The subscripts o2 represent the 
properties of the liquid oxidizer at the reference condition. The 
units used for the above expression (eq. 3.41) are defined in the list 
of symbols. 

The injection of the oxidizer, through the injectors which are 
rifled on the inside, produces a swirling effect which enhances the 
mixing of the oxidizer with the fuel. Therefore, eq. (3.41) is 
multiplied by 0.448 as recommended by Szuch. 
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The gas-phase mixing time delay ( ) is obtained from published 

experimental data correlated by: 


o = G 

m 


io'\i+i -u, 

c g 50 I 


(3.42) 


where 1 and 1 are the chamber and nozzle lengths. Since there is 
c n 

no nozzle in the preburner, 1^ is neglected resulting in a conserva- 
tive approach in computing o^. The exit gas velocity v^ is computed 
from: 

V = M VvR f (3.43) 

c c ^ g c 

and G is the functional relation plotted in Fig. 3.4 below. 



Figure 3.4 Gas-phase mixing delay based on experimentally 
observed chugging frequencies. 

Source: Szuch J. R., "Digital Computer Program for Analysis of 

Chugging Instabilities," NASA TN--D--7026, Dec. 1970. 
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The imaginary feed system impedances, and I^, must be 
computed, as well as the comparison of the critical real values, R^' 
and R^', with the operating point values. The computer program was 
utilized to calculate and at the operating point, solving 
equations (3.29) and (3.31) to obtain R^' and R^' at the specified 
frequency of interest. This calculation involves: 

1. breaking the feed system into elements, each having an 
impedance with real and imaginary parts a and P respectively. 

2. manipulating these elements into series and parallel combina- 
tions , and 

3. reducing the combinations to a single feed system impedance. 

The series combination of Z. and Z. (i.e. Z, =Z. +Z.) produces: 

1 J k 1 j ^ 

= a + a. (3.44) 

K I j 

= + (3.45) 

while the parallel combination of Z. +Z. {i.e. Z, = Z.Z./(Z.+Z.)} 

^ 1 J k 1 J ^ 1 

• produces : 


a.(a^ + P“) + a.(a^ + p“) 

Q _ J I J J 

(u . + a )^ + (P . + P 

‘ J .J 

P.(a^+ p^) + p.(a^ + P^) 
P = ■’ ^ I J J 

^ (a. + a.)2 + (P +p.)2 

I J 


(3.46) 


(3.47) 
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Using a lumped parameter approach, the feed system can be divided 
into the following elements: 

1. lines having only inertial pressure drop. The flow 
impedance therefore consists of a = 0 , and P = caZ/A , where 1 is the 
line length and A the cross-sectional area of the line. This is 
analogous to an electrical inductance 

2. orifices having only frictional pressure drop. The flow 
impedance therefore consists of 2AP/m^ , and P = 0 , where 

and m are the steady-state pressure drop across the injector element 
and mass flow rate respectively. This is analogous to an electrical 
resistance 

3. Storage volumes have an impedance consisting of a = 0 , and 

P = — B/o)pV > where B is the bulk modulus and V is the manifold 
S S 

volume. This is analogous to an electrical capacitance 

Fig. 3.5 shows the impedance representation for the oxidizer system 
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Fig. 3.5 Impedance representation of the oxidizer system. 
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The above system is collapsed to the following equation utilizing 
the series and parallel combinations. 


Z Z 

Z (s) = a + + 

“ Z +Z 

s man. 


(3.48) 


where a is the injector resistance ( ), b is the injector tube 

inductance (oLLA , and Z and Z are the suction line and manifold 

s man 

voltune impedance respectively. 

The real and imaginary parts of the oxidizer feed system 
impedance are computed according to the following: 


R = 2Ap/m 

O ^ O 


(3.49) 


and 


/ = cob + 

O 


((3-) + p if ) 

t s man 

fP^ + p2 ) 

* man 


(3.50) 


The imaginary part will be used to solve the characteristic 
equation and the real part will be compared to the critical value 
obtained from the solution of the characteristic equation. 

The same procedure holds for the fuel impedance Z^(s), although 
some simplifying assumptions are required to handle the system 
conveniently. Asscuning steady-state conditions upstream, the fuel 
injector may be considered to be fed by a constant fuel flow rate. 
Since the fuel is completely vaporized, the isothermal compressibility 
must be considered in computing Z^(s). With the above assumptions 
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the following equation is utilized to compute the fuel flow rate 
through the injector element. 



(3.51) 


where is the flow discharge coefficient. The flow discharge 
coefficient for the SSME fuel prebumer injector is calculated from 
conditions at rated power level and foiind to be 0.52. 

The steady-state fuel injection velocity v^, which is required 
for the computation of the drop size, is computed from: 



g f fi 
P A, 

Cl f 


(3.52) 


Linearizing and taking the Laplace transformation of eq. (3.51) 
results in the following equation. 


fh^. (s) = 


1 

P^fs) - 

P . 

Cl f 

P], - P^ 

Pl - P“, 

fl Cl 


A Cl 


P is) 

Cl 


(3.53) 


The perturbations in the injection manifold pressure may be expressed 
as : 


P,,(.) = 


R T 

g 

V 


0>0 

rh (s) 
A 

f S 


(3.54) 
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Combination of 


where is the fuel injector manifold volume, 
equations (3.52) and (3.53) and letting s = jco , results in: 


where 


(jw) = 


P. (s) - P is) 

ft Cl 

(s) 




pI - 

fl Cl 

P 

Cl f 


f 0)'P . 


R T 

g 


(3.55) 


(3.56) 


(3.57) 


The imaginary part, is used in the solution of the 
characteristic equation and the real part, R^, will be compared to the 
critical value obtained from the solution of the characteristic 
equation. 

The critical values of the pressure upstream of the injector 
element tube and the corresponding ratios of pressure drop to chamber 
pressure are computed from: 


AP 


R 'm 

o o 

2P 


IP. 


P -P 

fl Cl 


(3.58) 


(3.59) 
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Calculation of Stability Boundary 

A computer program capable of generating stability limits over a 
range of operating conditions has been developed and applied to the 
fuel preburner of the SSME. The program, drawing extensively from an 
earlier program by Szuch [12], is written in Fortran 77 and implemen- 
ted on the VAX 11/780 at the University of Tennessee Computing Center. 
The program solves the characteristic equation, utilizing the closed- 
form quadratic formula, for various sensitivity studies. This 
section will assist the reader in obtaining an overview of the compu- 
tation of the stability boundary. 

A range of possible chugging frequencies must be specified to 
initiate the program. It was found in the literature that the 
frequency range for chugging instabilities is between 75 to 200 Hz. 
Frequencies higher than 200 Hz will be neglected. At the selected 
operating point arid frequency of interest, the following computations 
are made by the computer program: 

1. calculate the steady-state value of imaginary fuel impedance 
from eq. (3.A7) 

2. calculate the real and imaginary oxidizer impedance and 
from equations (3.49 and 3.50) 

3. calculate the remaining variables in equations (3.29 and 

3.31) 

4. solve eq. (3.29) for 5^'* Negative values of R^' are 
neglected 

5. using the result from step (4), calculate R^' from eq. 
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(3.31). Negative values of R^' are also neglected 

6. convert results of R ' and R^' to critical values of 

of 

injection pressure drop from equations (3.58 and 3.59) 

7. store results from step (6) including the frequency and 
repeat steps (1) to (6) 

8. write out results from step (7) when frequency is out of . 
specified range 

A sorting routine is incorporated into the program to write out 
the critical pressure drop ratios ( ^ ^^fi^^c^ 

order of increasing ( AP^^/P^). 

The modified program is capable of producing sensitivity studies 
for variations in chamber pressure, fuel flow rate, oxidizer flow 
rate, fuel and oxidizer temperatures, and its affects on the 
instabilities are studied. The effects of these parameters on 
instability are presented in the following chapter. 
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CHAPTER IV 


MODEL APPLICATIONS AND RESULTS 


Verification of Nximerical Analysis 

The computer program, presented in Appendix B, was verified by- 
comparing the results obtained from the niomerical procedure to 
experimental results provided by NASA and Rocketdyne. The required 
program input values for full power level are presented in Table 4.1. 

The model was tested with values at full power level (FPL) and predicted 
stable operation in the fuel preburner of the SSME (Fig 4.1). This 
prediction conforms with experimental results since it was noted by 
NASA, during test stand firings and flight conditions, that the SSME is 
stable at steady-state conditions. 

Table 4.1 Input values at FPL obtained from Rocketdjme 
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FUEL INJ PRESS DROP TO CHAM PRESS 



Figure 4.1 Stable operation at full power level (FPL), 
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The analysis of the SSME fuel preburner chug during shutdown, 
presented in this thesis, is in agreement with the experimental 
results supplied. The frequency for the fuel preburner chug (Fig. 
4.2) varies from 115 Hz to 145 Hz, as determined by counting the 
peaks over a selected time segment. The mean frequency of the chug 
over the entire time interval is 125 Hz. The chug frequency 
predicted by therJiodel differed by 10 percent from the mean value 
obtained experimentally. This disagreement was due to the conser- 
vative linearized approach taken in the analysis. The predicted 
frequencies are in good agreement at the two end regions of the chug, 
since the low amplitude chug is more linear at the start and finish. 
The frequency at the two end points is approximately 140 Hz. Figures 
4.2 and 4.3, which are experimental results performed by NASA, show 
the chug frequency trace at the mean chugging pressure and the chamber 
pressure variation in the fuel preburner during shutdown respectively. 
These results, in conjunction with Figures 4.4 and 4.5, show the 
comparison between experimental and analytical predictions. 

In order to isolate the chug during experimental runs, the analog 
data were filtered to pass only frequencies on the range 70 < f < 200 
Hz and scanned visually for an amplitude surge representing the chug. 
When the surge was noticed the data were electronically plotted with 
pressure as a function of time. The fuel preburner chug trace 
presented in Fig. 4.2 shows that at a mean chamber pressure of 
4.4815E6 Pa (650 psia) the chug starts at approximately 2.3 
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Figure 4. 


Fuel Preburner chamber pressure variation (Coutesv of 
NASA Huntsville). 
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Figure 4.4 Instability of preburner. 

P = 4.4815E6- Pa, T = T. = 120 K 
c of 
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Figure 4.5 Stability of preburner with high oxidizer and fuel 

temperatures at maximxim allowable pressure . 

P = 3.7902E6 Pa, T = T. == 120 K 
c of 
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seconds after the shutdown command is issued and continues for about 


0.5 seconds before stabilizing out. 

The amplitude of the chug shown increases and then decreases with 
a pause at about 2.6 seconds before increasing again. This pheno- 
menon experienced by the chug is still unclear and further investiga- 
tion is needed in this area on amplitude prediction. 

Since the chug was noted experimentally to stable out at 
approximately 3.0 seconds after the shutdown command is given. Fig. 

4.3 shows the fuel preburner chamber pressure subsequently decreasing 
from 4.4815E6 Pa (650 psia) to about 6.8945E5 Pa (100 psia) after 3.9 
seconds. Therefore, it may be inferred that stable operation in the 
fuel preburner is prevalent in the region where the chamber pressure 
is below the mean chugging pressure of 4.4815E6 Pa (650 psia). 

The results obtained from the analytical solution, plotted in 
Figures 4.4 and 4.5, are in good agreement with the experimental 
results discussed above. The model was tested with a set of varying 
chamber pressures from 3.4474E6 Pa (500 psia) to 5.8605E6 Pa (850 
psia). The other required input conditions during the shutdown phase 
of the SSME is presented in Table 4.2. 

Figures 4.4 and 4.5 represent the nominal operating condition 
together with the critical pressure drop values ( 

( frequency. The only boundary point of interest 

is the critical value of the fuel injector pressure drop ( AP^^/P^). 

At pressures greater than or equal to 4.4815E6 Pa (650 psia), shown in 
Fig 4.4, the calculated operating point ~ 1*98 signifying 
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Table 4.2 Input values at shutdown condition 


Parameters 

Values 

P 

c 


T 

o 

120.0 K (216 R) 

Tf 

120.0 K (216 R) 

m 

o 




T 

1065.0 K (1917 R) 

c 


instability with a frequency of approximately 137 Hz. Fig 4.5 shows 
the marginal stability with an operating point, AP^^/P^ of 2.49, in 
the fuel preburner at a chamber pressure of 3.7902E6 Pa (550 psia). 
Because the operating point is just below the minimum point of the 
stability boundary where A P^^^/P^ is 3.27, stable operation is 
indicated. Therefore, it can be concluded from these results that 
the model predictions are in agreement with experimental techniques. 

Results of Sensitivity Study 

The onset of the chugs due to chamber pressure variations were 
analyzed utilizing the computer program. The chamber pressure was 
varied from 3.4474E6 Pa (500 psia) to 5.8605E6 Pa (850 psia) as input 
into the program. Since the chugging frequency ranges between 75 to 
200 Hz, the stability boundary computation is performed with the lower 
frequency bound and is incremented by 0.1 Hz to calculate the required 
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pressure drops. A stability boundary was therefore generated by 
computing the fuel injector pressure drop and the oxidizer injector 
pressure drop at the critical frequency of interest. The stability 
of the nominal operating point can thus be determined. 

The model was used to generate stability curves for several 
parametric studies. These studies include chamber pressure varia- 
tions, oxidizer and fuel temperature variations and oxidizer and fuel 
flow rate variations. The above parametric studies were also 

performed utilizing two different bulk moduli. The bulk moduli used 
in the model were for gaseous helium and liquid oxygen. Since helium 
is an inert gas the compressibility of helium may provide the. 
necessary softness for the chug. Unless otherwise stated, the value 
of the helium bulk modulus was used in the model. These parametric 
studies were performed with a single variable parameter while others 
were held constant. 

As indicated in Figures 4.4 and 4.5, instability exists at 
pressures greater than 3.7902E6 Pa (550 psia) while stable operation 
is prevalent at pressures below 3.7902E6 Pa (550 psia). At pressures 
greater than 4.4815E6 Pa (650 psia) instability exists, however, there 
is an upper bound to the instability since the preburner is stable at 
full power level. 

When low oxidizer and fuel temperatures (T^ = 40 K, T^ = 40 K 
(72 R)) were utilized as inputs, with varying chamber pressure the 
system was inevitably unstable even at low pressures with a frequency 
of 115 Hz and 121 Hz, shown in Figures 4.6 and 4.7 respectively. The 
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Figure A. 6 Unstable operation with low fuel and oxidizer 

temperatures at low chamber pressure. 
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Figure A. 7 Unstable operation in preburner with low oxidizer 

and fuel temperatures at mean chugging pressure . 

T = = 40 K, P = A.A815E6 Pa 

of c 
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ratio of fuel injector pressure drop to chamber pressure ( 
of the operating point was never below the minimum point of the 
stability boundary. 

The effect of fuel and oxidizer temperatures on the stability of 
the system was also considered. At a chamber pressure of 4.A815E6 Pa 
(650 psia), which is the mean chugging pressure, the oxidizer 
temperature was subjected to a variation ranging from 40 K to 120 K 
while the fuel temperature was kept relatively low at 40 K. The 
variation of varying oxidizer temperature with low fuel temperature 
(Fig. 4.8) did not stabilize the fuel preburner system, however, it 
was noted that operation at high fuel and oxidizer temperatures (120 K 
and 100 K respectively) stabilized the system shown in Fig. 4,9. 

High oxidizer and fuel temperatures cause partial vaporization of the 
liquid oxidizer; hence smaller droplet radius. Vaporization and 
mixing times are thus decreased with the smaller droplet radius. The 
stability of the system, probably due to the proper mixing of the 
reactants, could be achieved by increasing propellant temperatures 
even at the mean pressure where chugging is prevalent. 

Since it is undesirable to work with high oxidizer temperatures, 
several test cases were run to determine if oxidizer temperatures 
could be reduced and yet maintain stability. It was noted that 
stability was achieved with higher fuel temperatures while other 
parameters remained constant- This result presented in Fig. 4.10 
shows the fuel temperature at 150 K while oxidizer temperature was 
lowered to 86 K. 
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Figure 4.8 


Unstable operation in prebumer with low fuel and 
high oxidizer temperature at mean chugging pressure, 
T = 40 K, T^- = 120 K, P = 4.4815E6 Pa 
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Figure A. 9 


Stable operation in preburner with high fuel and 
high oxidizer temperatures at mean chugging pressure. 


T = 100 K, T. = 120 K, P = A.A815E6 Pa 
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Figure 4.10 Stable operation in preburner with increasing fuel 
temperature and lower oxidizer temperature at mean 
chugging pressure. 

T = 85.7 K, T. = 150 K, P = 4.4815E6 Pa 
o I c 
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At the nominal operating point, reduction of fuel flow rate is 
unlikely to inhibit the chug, since chugging exists at relatively high 
fuel flow rate (equivalence ratio = 6.0). The mean droplet radius, 
being an inverse function of fuel flow rate, increases with fuel flow 
reduction. This causes an increase in vaporization and mixing time 
and hence increases the combustion instability. Figures 4.11 and 
4.12 show the operating point in the unstable region with fuel flow 
reduction at a frequency of 129.2 Hz and 136.3 Hz respectively. 

The influence of varying the oxidizer flow rate on the chug was 
studied. Although there is no oxidizer feed during shutdown, as the 
FPOV is already closed, the helium flow rate is considered as the 
oxidizer flow since it clears the residual oxidizer from the lines and 
manifold. This is possible because the residual oxidizer flow rate 
is equivalent to that of helivun using the appropriate density value. 
The model predicted that stability was achieved at low pressures with 
increased flow rates. 

The model was also tested with the bulk modulus of liquid oxygen 
at conditions where chugging exists. There were no stability bounda- 
ries generated for those conditions signifying stable operation. 

The frequency range of 75 to 200 Hz was tested with variations in the 
chamber pressure, oxidizer and fuel temperatures. The instability 
which was prevalent at previous operating conditions, with the bulk 
modulus of helium, was eliminated utilizing the bulk modulus of liquid 
oxygen. This shows that the chug depends not only on operating 
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Figure 4.11 Unstable operation in preburner with low fuel flow 

rate at mean chugging pressure. 

m. = 5 Kg/ sec, P = 4.4815E6 Pa 
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Figure A. 12 Unstable operation in preburner with low fuel flow 
rate at mean chugging pressure. 
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conditions but also on the helium purge conditions, specifically the 
compressibility of the helium. 

The following figures (Figures 4.13, 4.14 and 4.15) show the 
variation of chugging frequency with oxidizer temperature, fuel flow 
rate and chamber pressure respectively. Fig. 4.13 shows that stable 
operation was achieved with oxidizer temperatures between 80 to 115 K. 
Temperatures below 80 K resulted in unstable operation at a frequency 
of 139 Hz. -F^ig. 4.14 shows that the system was unstable with fuel 
flow variation at a frequency of 138 Hz. The two endpoints of the 
curve were arbitrary selected for this parametric study. Fig. 4.15 
shows that stable operation is permissible at low pressures. The 
chugging frequency increases somewhat linearly with increased, chamber 
pressure up to 5.759E6 (830 psia). 

The preceding parametric studies reflected several intriguing 
results. Although the analysis in this model was simplified, never- 
theless, it provided the basic understanding to the occurrence of the 
SSME fuel preburner chug and the parameters to be investigated further 
for chug elimination. 
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Figure 4.13 Variation of chugging frequency with oxidizer 
temperature . 
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Figure 4.14 Variation of chugging frequency with fuel flow 
rate. 
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Figure 4.15 


Variation of chugging frequency with chamber 
pressure. 
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CHAPTER V 


CONCLUSIONS AND RECOMMENDATIONS 


The objectives of this research were to identify the significant 
elements in triggering the onset of the chug and its' possible 
elimination. System variables such as chamber pressure, propellant 
temperatures and flow rates were varied and their effects -on the 
stability of the SSME fuel preburner chug analyzed. 

The first task was to review the literature in the area of low 
frequency instability as well as to solicit available computer models 
capable of analyzing the chug. The best available program was modi- 
fied to analyze the fuel prebumer chug. The model presented in this 
thesis is not restricted to SSME preburners but may also be applied 
to liquid propellant rocket engines. This model proved to be 
valuable in meeting the required objectives. 

The model in this thesis assumes quasi-steady state conditions 
and was linearized to simplify the non-linear governing differential 
equation. The method of perturbations was used to linearize the 
governing differential equation. Chemical kinetics in the combustion 
process were assimied to be infinitely fast and were therefore 
neglected in the analysis. 

Although the model has been simplified and linearized, it is 
nevertheless capable of chug predictions and is in excellent agreement 
with experimental results provided by Rocketdyne and NASA. The 
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predicted frequencies differing by less than ten percent of the 
measured value show that this model is effective in analyzing the SSME 
fuel preburner chug. 

The VAX 11/780 CPU time required for each sensitivity study was 
approximately 5 seconds. This extraordinarily short GPU time is 
beneficial as several different test cases could be performed with 
limited cost compared to experimental runs. 

Since the chug occurs during the heliiim purge of the oxidizer, 
the sensitivity studies were performed utilizing the bulk modulus of 
heliiam and repeated with the bulk modulus of the liquid oxidizer. 

The chugs did not exist when the relatively incompressible oxygen was 
considered. 

The results of these sensitivity studies are presented and 
discussed in detail in the previous chapter. Table 5.1 summarizes 
the results. Several intriguing results were revealed by the 
model. The onset of the chugs were linked to a couple of 
significant parameters such as fuel temperatures and fuel flow 
rates. The predicted results show that the fuel preburner 
instability was suppressed at higher fuel temperatures. This 
explains why chugging is reduced at high cut-off power; hence high 
fuel temperature. It was also noted that at high fuel 
temperatures the oxidizer temperature could be reduced. Other 
parameters of interest such as oxidizer flow rates and bulk modulus 
have shown to influence the instability to some extent. 



Table 5 . 1 Summary of parametric studies performed 


Variation of System 
Variables 

Comments 

P 

c 

At low fuel and oxidizer temperatures 
approximately 40 K (72R) the system was 
inherently unstable. At the other extreme 

with high fuel and oxidizer temperature 
= 120 K (216 R) stability is 
achieved with maximum pressure of 3.7902E6 
Pa (550 psia) 

T 

o 

At the mean chugging pressure of 4.4815E6 
Pa (650 psia) low oxidizer temperature and 
high fuel temperature proved unsuitable for 
operation. However, with increased 

T = 100 K (180 R) stability is possible. 
Further increases in fuel temperature to 
150 K (270 R) showed that stability was 
achieved with lower oxidizer temperatures. 

Tf 

At low fuel temperatures the chug was not 
inhibited at all with varying oxidizer 
temperatures even at low pressures where 
stability was prominent. 

m 

o 

High helixim flow rates (analogous to 
oxidizer flow) showed that chugging was 
eliminated. 

mf 

Variation in fuel flow rate did not 
eliminate the chug. It was noted that 

since the MFV is not closed until after 
chugging ends, fuel variation would 
not eliminate the chug. 

Bulk Modulus 

When liquid oxygen bulk modulus value 
was used the system was inherently stable 
since no stability boundary was generated. 
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The modeling effort should be continued and supported in an 
attempt to understand the non-linear chugs. This model, being 
linearized is limited in the capability of predicting amplitude 
changes of the chug. The challenge now is to incorporate the 
NASA/Rocketdyne transient model into the analysis. 

This research has laid the foundation for continuing study of the 
SSME shutdown chug. Although this model was developed to analyze the 
fuel preburner chug it may be applied to the oxidizer preburner or 
other liquid propellant rocket engine. This model will assist in the 
understanding of the chugging phenomenon as it was developed 
predominantly as a useful analysis tool for engine designs. 
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APPENDIX A 

DERIVATION OF EQUATIONS 
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In this appendix the reduction process involved in acquiring eq. 
(3.5) from eq. (3.3) as well as the general reduction of the charac- 
teristic equation will be presented. The method of perturbations 
will be utilized to linearize any non-linear differential equations. 
Eq. (3.3) is reproduced below; 


6 

dt 





(3.3) 


Taking the perturbations of system variables and neglecting higher 
powers and products of perturbations results in the following 
equation . 



p + p = 


m^(t) 


m^(t) 






m (f) 

e 


th it) 

e 


(Al) 


Dividing both sides by the total propellant flow rate gives: 


a 




m m At) + m At) — rh (i) 
, p fh ob e fh ob e 

-!! + == — : + : 

0 ) 


(A2) 


which can be reduced easily to; 


B 

9 - 
- Bt 



m At) + rh At) 
;b ub 

m 

t 


(3.5) 
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Combination and reduction of characteristic equation 

The characteristic equation, derived in chapter three, is 
presented below as: 


P XK , P YK , 
_£ — ^ o ^ _£ — ^ r 




Z„u-) 


Z^(s) 


0 s + Cl 

8 


(3.24) 


Letting s = jco Z^(s) = R^' + Zj^(s) = R^' + and noting that 


■IX 


= cos(x) - isin(x), eq. (3.24) thus becomes: 


-Cl -Jco'9^ = 


P XKcosioi’i ) ~jP XK sin (co'i ) 

c 0 c a 

R'+Jl 

0 o 


P V'/f cos (co'tj — jP YKsin (ca'tJ 
_£ r c f_ 

r;+ji. 

t t 


(A3) 


Multiplying both sides of eq. (A3) by (R^' + jl^) and (R^' + jl^) and 
separating into the respective real, and imaginary parts the following 
equations are obtained: 
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Real Part 


R 'R.’Cl - u'O (/ R : 

or s o r 


R 'I.) - [ [.Cl = PXK 

Of « r 


R cosiur'i ) + sm (o)'c ) 

r » r o 


+ PYK 


R ' co.'iiiji’z.) + [ sin {(jo’i.) 
» for 


(A4) 


Imaginary Part 


Cl 

\r ’[. + r;i 

+ co’0 

[ [,-R 'R : 

\ = p xk\ 

7. cos (co't ) — R sim (w'z ) 


[or to 

1 

of or 

1 ^ 1 

r Of 0 


+ P YK 

c 


I cos (oo't: ,) — R ' sin (co't .) 
0 f o ' r 


(A5) 


Solving equations (AA and A5) for R^' ; 
Real Part 



P YK 

c 


R 'cos (co't J + ^ sin (co't J 

o ^ f o f 


+ IJ^ XK sin (co't ) - co'e R + ClI I, 

re o g o f of 


R 'Cl + co'0 I — P XK cos (co't ) 

0 g o c o 


(A6) 


Imaginary Part 



P YK 


R 'sin(o)’iJ — I cos (co't.) 
o tor. 


- [^ XKcos(ui\ ) + co'9 / /, + CIR I. 
r e o got or 


R ’(o'0 -CU -P XKsiniw'i ) 

o g o c o 


(A7) 
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Equating equations (A6 and A7) and solving for R^' yields the 
characteristic equation in quadratic form. 



o 



P YK 

c 


Cl- + co'“0^ 

g 


+ 


R 


2XP 


Cl cos ( o)'e ) + CO '6 sin (co't ) 
o g o 


P XK sin (co'c 

c o 


— co’r^) 


I^F ~ I P XKcosioi’i - coV) - 

0 0 c o f 



P YK 

c 


(C1- + co’‘^0-) 
g 


PP X^K 

f C 

Y 


2X1 [, 

O ! 

~ 


Cl sin (co 't ) — 0 ) ' 0 cos (to 't ) 

0 3 0 


= 0 


(3.29) 


where 


F = CO '9 cos((o ’ c.) — Cl sin (to ' cj 

g r r 
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APPENDIX B 

LISTING OP PROGRAM CHUGTEST 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


iHHHnHnHririhnfinr 


# 

ir^ 

# 

ir4nrinriririnrjririnrin,-iririnritir-;nrmnrirViririr4nni-ir'iriririr1rir1ririi-!rifTrinrinrinf7i 


PROGRAM CHUGTEST 
BY KAIR-CHUAK LIM 


# 

# 


ir^ 
# 
# 
il: 


THE PROGRAM BELOW WAS DEVELOPED TO ANALYSE THE 
COMBUSTION INSTABILITY (CHUGGING) IN THE FUEL PREBURNER 
OF THE SSME. 


# 
# 
// 
# 
# 


DIMENSION UPL(IOOO), ACRL(IOOO), UPH(IOO) , ACRH(IOG) , FL(IOOO), 

* FH(IOO) 

DIMENSION ACCL(IO) ,UPPL(10) , FFL(IO), ACCH(IO), UPPH(IO), 

* FFH(IO) 

REAL NUMER 

INTEGER ENO, ENF, PARAM 

OPEN (UNIT=20 , FILE= ' OUTCHUG . DAT ' , STATUS= ' OLD ’ ) 

OPEN (UNIT=21,FILE=' PL0T21.DAT' , STATUS= ' OLD ' ) 

OPEN (UNIT=22,FILE=' PLOT22.DAT' , STATUS= ' OLD ' ) 

OPEN (UNIT=23 , FILE= ' PL0T23 . DAT ' , STATUS= ' OLD ’ ) 

OPEN (UNIT=24,FILE=' PL0T24.DAT’ , STATUS= ' OLD ’ ) 

OPEN (UNIT=25,FILE=’PLOT25.DAT' , STATUS= ' OLD ' ) 

OPEN (UNIT=26,FILE=' PLOT26.DAT' , STATUS= ' OLD ’ ) 

OPEN (UNIT=27 ,FILE=' PLOT27 .DAT' , STATUS= ' OLD ' ) 

OPEN (UNIT=28,FILE=' PL0T28.DAT' , STATUS= ' OLD ’ ) 

IOUT=20 
10=21 

Cl=-0.1072 
PI=4.*ATAN(1.) 

1 FORMAT (//) 

C THE FOLLOWING DATA IS FOR 
RF=4126. 

RH0F=20 . 

C THE FOLLOWING DATA IS FOR 
TCR=154.5 
WT0=32. 

HV=2.13E5 
RH0=1152. 

ST=1.19E-2 
VIS=3.0348E-4 

RSTVC=((1143./RH0)*(VIS/1.9E-4)*(ST/1.33E-2))**.25 
WTT=(WTO/100. )**.35 ! MOL'.- WT COEFF USED IN EL50 


HYDROGEN FUEL 

! FUEL GAS CONSTANT (J/kg C) 
! FUEL DENSITY (kg/m3) 

CYROGENIC OXYGEN OXIDIZER 
! CRITICAL TEMP (K) 

! MOL. WT (kg/kg mole) 

! HT. OF VAPORIZATION (J/kg) 
! DENSITY (kg/m3) 

! SURF. TENSION (N/m) 

! VISCOSITY (N s/m2) 
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HVT=(HV/3.26E5)**.8 
PHI=6 . 

SWC=0 . 448 


! HT. OF VAP COEFF USED IN EL50 
! EQUIVALENCE RATIO 
! SVRILING FACTOR 


C THE TYPE OF PARAMETRIC STUDY IS SELECTED AT THIS POINT 
15 WRITE (5, 2) 

2 FORMAT (IH ,'TYPE OF PARAMETRIC STUDY? 3x 1 CHAMBER 

@ PRESSURE ' , / , 3x , ' 2 OXIDIZER FLOW RATE ' . / , 3x , 
ta '3 FUT:L FLOW RATE' ,/,3x, '4 OXIDIZER TEMPERATURE , 

@ 3x , ' 5 FUEL TEMPERATURE ',///, 3x , ' ENTER SELECTION ' ) 

READ (5, *) PARAM 


C THE FOLLOWING IS COMBUSTION CHAMBER GEOMETRY DATA OBTAINED 


C FROM NASA AND ROCKETDYNE 
CDIA=0 . 265 
ELC=0 . 108 
AC=PI*CDIA**2/4 . 
VC=AC*ELC 

C THE FOLLOWING IS SPECIFIED 
. D0=2.8956E-3 
ENE=264 . . 

ENO=ENE 

ENF=ENE 

ELO=0.0511 


! CHAMBER DIAMETER (m) 

! CYL. CHAMBER LENGTH (m) 

! COMB CHAMB. X-AREA (m2) 

! CHAMBER VOL (m3) 

INJECTOR GEOMETRY 
! OX INJECTOR ORIF. DIA (m) 
! NO. OF ELEMENTS 
! NO OF OX ELEMENTS 
! NO OF FUEL ELEMENTS 
! OX ELEMENT LENGTH 


C FOR SINGLE ORIFICE INJECTOR, VOM IS OXIDIZER INJECTOR MANIFOLD 
C VOLUME. VOF IS FUEL MANIFOLD VOLUME 


CD0=.5 

CDF=.5 

SWC=.448 

V0F=6.55E-5 

ELLS=0.106-' 

DSU=0 . 0445 

V0M=9.34E-4 

APF0=1.329E-5 


! FLOW COEFF OX INJ ELEMENT 
! FLOW COEFF FUEL INJ ELEMENT 
! CORRECTION FOR SWIRLER 
! FUEL INJ MANIFOLD VOL (m3) 

! SUCTION LINE LENGTH (m) 

! SUCTION LINE DIA (m) 

! OX MANIFOLD VOL (m3) 

! FUEL ANNULUS (m2) 


C AT THIS POINT A DECISION IS MADE WHETHER TO GENERATE BOUNDARIES 
C AT VARIOUS CHAMBER PRESSURES, USING OPERATING POINT DELAYS OR 
C TO GENERATE BOUNDARIES USING DELAY VALUES SENSITIVE TO FUEL 
C ANNULUS AREA JJ=0 DENOTES FORMER AND JJ=1 DENOTES L.\TTER 


J=1 

JJ=0 

MN=1 

LI=0 

IL=0 

IF (JJ .EQ. 0) MN=18 


TYPE OF STUDY 
LOWER THROTTLING LIMIT 
LOW FREQ SOLN INDEX 
HIGH FREQ SOLN INDEX 
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C THIS SECTION DIVIDES THE RANGE OF THE PARAMETRIC STUDY INTO 
C EIGHT EQUAL SEGMENTS 

IF (PARAM .EQ. 1) THEN 
WRITE (5, 3) 

3 FORMAT (IH , 'WHAT IS THE MINIMUM AND MAXIMUM PRESSURE RANGE 
@ FOR STUDY?') 

READ (5, *) PCI, PC2 

COUNT= ( PC 2 - PC 1 ) / (FLO AT ( MN ) - 1 . ) 

WRITE (5, 4) 

4 FORMAT (IH ,' INPUT OXIDIZER AND FUEL FLOW RATE AND 
@ TEMPERATURE RESPECTIVELY') 

READ (5, *)WO, WF, TO, TF 
GO TO 13 

ELSE IF (PARAM .EQ. 2) THEN 
WRITE (5, 5) 

5 FORMAT (IH ,'WHAT IS THE MINIMUM AND MAXIMUM OXIDIZER FLOW 
@ RATE RANGE FOR STUDY') 

READ (5, *) WOl, W02 

COUNT= (W02 -WOl)/ (FLOAT (MN ) - 1 . ) 

WRITE (5,6) 

6 FORMAT (IH ,' INPUT CHAMBER PRESSURE, OXIDIZER AND FUEL 
@ TEMPERATURE FUEL FLOW RATE RESPECTIVELY') 

READ (5, *) PC, TO, TF, WF 
GO TO 13 

ELSE IF (PARAM .EQ. 3) THEN 
WRITE (5, 7) 

7 FORMAT (IH , 'WHAT IS THE MINIMUM AND MAXIMUM FUEL FLOW RATE 
@ RANGE FOR STUDY') 

READ (5, *) WFl, WF2 
COUNT=(WF2-WFl)/ (FLOAT(MN)-l. ) 

WRITE (5, 8) 

8 FORMAT (IH 'INPUT CHAMBER PRESSURE, OXIDIZER AND FUEL 
@ TEMPERATURE AND OXIDIZER FLOW RATE') 

READ (5, *) PC, TO, TF, WO 
GO TO 13 

ELSE IF (PARAM .EQ. 4) THEN 
WRITE (5, 9)’ 

9 FORMAT (IH , 'WHAT IS THE MINIMUM AND MAXIMUM OXIDIZER 
@ TEMPERATURE RANGE FOR STUDY') 

READ (5, *) TOl, T02 

COUNT= (T02 -TO 1 ) / (FLOAT (MN ) - 1 . ) 

WRITE (5, 10) 

10 FORMAT (IH ,' INPUT CHAMBER PRESSURE, OXIDIZER AND FUEL FLOW 
(a RATE AND FUEL TEMPERATURE') 

READ (5, *)PC, WO, WF, TF 
GO TO 13 

ELSE IF (PARAM .EQ. 5). THEN 
WRITE (5, 11) 

11 FORMAT (IH ,'WHAT IS THE MINIMUM AND MAXIMUM FUEL 
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12 


13 


@ TEMPERATURE RANGE FOR STUDY') 

READ (5, *)TF1, TF2 

COUNT= (TF2 -TF 1 ) / (FLOAT (MN ) - 1 . ) 

WRITE (5, 12) 

FORMAT (IH INPUT CHAMBER PRESSURE, OXIDIZER AND FUEL 
@ FLOW RATE AND OXIDIZER TEMPERATURE') 

READ (5, *)PC, WO^ WF, TO 
END IF 

DO 90 1=1, MN 
INDEX=I-1 

IF (PARAM .EQ. l)PC=PCl+COUNT’'^FLOAT( INDEX) 

IF (PARAM .EQ. 2)W0=W01+C0UNT>^FL0AT(INDEX) 

IF (PARAM .EQ. 3)WF=WF1+C0UNT*FL0AT(INDEX) 

IF (PARAM .EQ. 4)T0=T01+C0UNT*FL0AT(INDEX) 

IF (PARAM .EQ. 5)TF=TFl+COUN'T‘'^FLOAT(INDEX) 


OF=WO/WF 

BO=PC 

K=0 

APF=APFO 

AOV=PI*0.0279^' 


^2/4. 


MIXTURE RATIO 
BULK MODULUS (Pa) 
VARIABLE FUEL AREA INDEX 
FUEL ANNULUS AREA (m2) 
CTRL VAL AREA (m2) 


C COMPUTE SLOPE OF X AND Y FROM EQUATIONS ' 

SLOPE=l . 3644*PHI**3-38 . 7501*PHI**2+376 . 18*PHI-1340 . 08 


C THE SPECIFIC HEAT RATIO IS USED TO CALCULATE CHAMBER PRESSURE 
C AT THE INJECTOR FACE 

TC=1065. ! COMBUSTION TEMP 

RC=122.3 ! COMB PROD GAS CONST (J/Kg K) 

GAMMA=0. ! SPECIFIC HEAT RATIO OF PRODUCTS 

CALL GAM (TC, GAMMA) 

ACOUS=SQRT(GAMMA*RC*TC) ! ACOUSTIC VEL EXIT 

vrroT= WO + wf 

RHOH2=PC*l .E-6*2 . / (82. 06E-3*TC*101325*1 . 03) 

RH0H20= ( 1 . / 1 . 834) / (2 . 2*2 . 54E-2**3*12 . **3) 

RH01=(5./20.)*RHOH2 + (1./20. )*RH0H20 
AE=4.92E-2 

EMC=WT0T/(AC0US*RH01*AE) !MACH NO AT EXIT 

VELC=EMC*ACOUS ^ ! GAS VEL EXIT 

EXP 1=GAMMA/ (GAMMA - 1.) 

PC I=PC* ( 1 . +GAMMA*EMC**2 )/(!.+ (GAMMA- 1 . ) / 2 . *EMC**2 ) 

@ **EXP1 ! INJ. PRESS 

C COMPUTE VELOCITIES AND PRESSURE DROPS. 

14 VF=WF/ENF/PCI*RF*TF/APF ! FUEL INJECTION VEL 

PF=PCI*SQRT(1.+VF**2/CDF**2/RF/TF) ! FUEL INJ MANIFOLD PRESS 
V0=W0/EN0/RH0/PI/D0**2*4. ! INJ VEL OF DROPLET 
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c 


EMR=OF*VO/VF 

DPJ=(WO/ENO/CDO/PI/DO* 

DPS= ( WO / AOV ) **2 / RHO/ 2 . 

DPP=DPJ+DPS 

PO=PCI+DPP 


! INJ MOMENTUM RATIO 
2*4.)**2/rhO/2. ! OX INJ-ORI PRESS 

DROP 

! PRESS DROP SUCTION LINE 
! PRESS DROP VALVE TO INJ FACE 
! OX MAN PRESSURE 


C THE FOLLOWING SECTION COMPUTES DROP SIZE, BURNING LENGTH, 
C DELAY TIMES. 


17 RMCTN=0.118*D0*SWC*RSTVC*SQRT(EMR) ! MEAN OX DROP RAD 

TT=(1.-T0/TCR)**.4 ! TEMPERATURE COEFF. 

ELS 0=0.069 (RMCTN/ 7 . 62E -5 )** 1 . 45* (VO/ 30 . 5 ) ** . 7 5* 

@ WTT*HVT/(PC/2.07E6)**.66 ! LENGTH FOR 50% PROPELLANT 

C VAPROIZED (id) 

TAUV=EL50/V0 ! OX VAP TIME DELAY (sec) 

WT0T=WF+W0 ! TOTAL PROPELLANT FLOW 

C RATE (kg/ sec) 

C COMPUTE GAS RESIDENCE TIME, GAS PHASE DELAY COEFF AND OPERATING 
c POINT PRESSURE DROP RATIOS. 

20 THETAG=RHOF*VC/\7TOT ! GAS RESIDENCE TIME (sec) 

DPO=(PO-PCI)/PC ! OX-INJ ORI PRESS DROP TO CHAMB PRESS 

DPF=(PF-PCI)/PC ! FUEL-INJ PRESS DROP TO CHAMB PRESS 

C0EFF=1.E3*(ELC-EL50)/VELC ! GAS PHASE DELAY COEFF 

TAUM=TM (COEFF) ! GAS PHASE MIXING AND RX TIME DELAY 

TAUT=TAUM+TAUV ! TOTAL OX TIME DELAY 

C THE FOLLOWING SECTION COMPUTES LINEARIZED RESISTANCES, 

C INDUCTANCES, CAPACITANCES AND GAINS 

XG=((1.+(0.5*(1.+PHI)/TC)*SLOPE)/WTOT)*(PCI/PC) 

FG= ( 1 . - ( 0 . 5* ( PHI* ( 1 . +PHI ) /TC ) *SLOPE ) / WTOT ) * ( PC I / PC ) 
RESF=(PF**2-PCI**2)/(PCI*WF) ! LINEARIZED FUEL 

ELJ=EL0*4./EN0/PI/D0**2 ! OX-INJ ELEM INDUCTANCE (1/ra) 
ELSU=ELLS*4./PI/DSU^~’^2 ! SUCTION LINE INDUCTANCE (1/m) 
COM=-BO/(RHO*VOM) ! OX INJ MAN. CAP. (1/sec m2) 

COF=VOF/RF/TF ! FUEL INJ MAN. CAP (1/sec m2) 

C PARAMETERS OF INTEREST ARE WRITTEN OUT 

WRITE (lOUT, 21) 

21 FORMAT (IH ,4X , ' PC ' , lOX, ' PO' , lOX, ' PF' , 7X, ’MO’ , 6X, ’ MF ’ , 3X, 

@ 'FTEMP' ,2X, 'OXTEMP' ,3X, 'TETAG' ,5X, 'TAUM' ,5X, 'RADIUS' ,6X, 

@ 'DPO' ,9X, 'DPF' ) 

WRITE (lOUT, 22) 

22 FORMAT (IH ,3X, ' (Pa) ' ,8X, ' (Pa) ' ,8X, ' (Pa) ’ ,4X, ' (Kg/s) ' ,2x, 

@ ' (Kg/s) ’ ,2x, ' (K) ' ,4x, ’ (K) ' ,5x, ' (sec) ' ,5x, ' (sec) ' ,6x, ' (m) ' , 

@ 7x,'(Pa)',8x,'(Pa)’,/) 
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WRITE (lOUT, 23) PC, PO, PF, WO, WF, TF, TO, THETAG, TAUM, 

@ RMCTN, DPO, DPF 

23 FORMAT (IH ,E10 . 4 , 2X,E10 .4 , 2X,E10 .4 , 2X,F5 . 2 , 2X,F5 . 2 , 2X,F5 . 1 , 
@ 2X,F5.1,2X,ES.2,2X,E8.2,2X,E8.2,2X,E10.4,2X,E10.4,//) 

C THE FREQUENCY BANDS AND INCREMENT ARE SELECTED 

C AND THE SOLUTION OF THE CHARACTERISTIC EQUATION IS INITIATED 

FREQ=74.9 

JN=0 

N=0 

NN=0 

27 FREQ=FREQ + .1 

IF ((FREQ .GT. 200.) .AND. (FREQ .LT. 300.)) GO TO 28 
GO TO 29 

28 IF (N .EQ. 0) JN=1 

N1=N 

FREQ=300 . 

JM=0 

N=0 

NN=1 

GO TO 30 

29 IF (FREQ .GT. 300.) GO TO 50 

30 W=6.2832*FREQ ! ANGULAR FREQUENCY 

11=1 ! MULTIPLE OX SOLN INDEX 

RX=0. ! CRIT. REAL PART OF OX IMPEDANCE 

RX1=0 . 

RX2=0 . 

RRF=0. ! CRIT. REAL PART FUEL IMPEDANCE 
RRF1=0. 

RRF2=0 . 

DP0C2=0. ! RATIO OF OX-INJ PRESS DROP TO CHAMB PRESS 

DP0C5=0 . 

DP0C8=0. 

DPFC=0. ! CRIT VAL OF FUEL-INJ PRESS DROP TO CHAMB PRESS 
DPF 1=0. 

DPF2=0 . 

C TERMS IN THE CHARACTERISTIC EQUATION ARE EVALUATED AT THE 
C SPECIFIED FREQUENCY 

EMAGF=-1./W/C0F*PF/PCI ! IMAG PART OF FUEL IMPEDENCE 

THETAT=TAUT*W 

THETAM=TAUM*W 

THETAV=TAUV*W 

THETAL=ELJ*W 

ALPHAG=THETAG*W 

Al=COS(THETAM) 

B1=SIN(THETAM) 

A2=C0S(THETAV) 
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B2=SIN(THETAV) 
A3=C0S (THETAT) 
B3=SIN(THETAT) 
F=ALPHAG*A1-C1*B1 


C THE SECOND AND FIRST ORDER COEFF OF THE CHARACTERISTIC 
C QUADRATIC EQUATION FOR CRITICAL OX REAL PART 
C ARE EVALUATED 

AG=F-EMAGF/(PC*FG)*(C1’'”'^2+ALPHAG**2) 

BG=-PG*XG*B2+2.*XG*EMAGF/FG*(C1*A3+ALPHAG*B3) 

C THE FOLLOWING SECTION COMPUTES REAL AND IMAGINARY PARTS OF THE 
C IMPEDANCES LOOKING UPSTREAM FOR A SINGLE ORIFICE CONFIG. 

AXA=2.*(PCI-P0)/W0 ! REAL PART OX-IMP (N s/ml kg) 

NUMER=((C0M/W)*(ELSU*W)**2 + (ELSU*W)*(C0M/W)**2) 
DENOM=(ELSU*W + C0M/W)**2 
BXA=NUMER/DENOM 

32 EMAGO=THETAL+BXA ! IMAG PART OF OX IMPEDANCE 

C THE ZEROTH ORDER COEFFICIENT IN THE QUADRATIC EQUATION FOR 
C CRITICAL OX REAL PART IS EVALUATED. THE FOLLOWING SECTION 

C SOLVES THE QUADRATIC AND RELATES SOLUTION TO THE PRESSURE 

C DROP RATIO. 

CG=EMAG0**2*F-EMAG0*PC*XG*A2-2.*XG*EMAGF*EMAG0/FG* 

@ (C 1*B3 -ALPHAG*A3 ) -EMAGF*EMAGO**2/ (PC*FG) * (C 1**2+ALPHAG**2 ) 
@ -XG**2*PC*EMAGF/FG 

C, THIS SECTION STARTS THE COMPUTATION OF THE CRITICAL 
C REAL VALLTIS OF THE OXIDIZER IMPEDANCE (Ro) 

IF (AG .NE. 0.) GO TO 33 
RX= -CG/BG 

RPC=RX-AXA ! CRIT REAL PART OF OX-INJ IMP 

IF (RPC .LE. 0.) GO TO 27 
DPOC2=RPC*WO/PC/2. - 
GO TO 38 

33 ALPHA=(BG**2-4.*AG*CG) ! SQRT FACT IN QUADRATIC EQ. 

IF (ALPHA) 27, 34, 35 

34 RX=-BG/2./AG 
RPC=RX-AXA 

IF (RPC .LE. 0.) GO TO 27 
DP0C2=RPC*W0/PC/2 . 

GO TO 38 

35 RXl=f-BG+SQRT(ALPHA))/2./AG ! CRIT VAL REAL PART OX-INJ 

RX2= ( - BG - SQRT (ALPHA ) ) / 2 . / AG 

RPC1=RX1-AXA 

RPC2=RX2-AXA 


104 



IF ((RPCl .LE. 0.) .AND. (RPC2 .LE. 0.)) GO TO 27 
IF (RPCl .LE. 0.) GO TO 36 
IF (RPC2 .LE. 0.) GO TO 37 
11=2 

DPOC5=RPCl*WO/PC/2. 

DP0C8=RPC2*W0/PC/2 . 

GO TO 39 

36 RPC=RPC2 
DP0C2=RPC*W0/PC/2. 

R.X=RX2 ■ ' 

GO TO 38 

37 RPC=RPC1 
DP0C2=RPC*W0/PC/2. 

RX=RX1 

GO TO 38 

C THE CRITICAL FUEL REAL PART IS DETERMINED USING OX SOLUTIONS. 

C THE RESULTS IS RELATED TO THE PRESSURE DROP 
C RATIO 

38 • RRF=(FG*RX*B1*PC-FG*EMAG0*A1*PC-XG*EMAGF*A3*PC+ALPHAG*EMAG0* 

@ EMAGF+EMAGF*RX*C 1 ) / (RX*ALPHAG-EMAGO*C 1 -XG*B3*PC ) 

IF (RRF) 27, 27, 43 

3 9 RRF 1= ( FG*RX 1*B 1*PC -FG*EMAGO*A 1*PC -XG*EMAGF*A3*PC+ALPHAG*EMAGO* 

@ EMAGF+EMAGF*RX1*C1) / (RX1*ALPHAG-EMAG0*C1-XG*B3*PC) 

RRF2= (FG*RX2*B 1*PC -FG*EMAG0*A1*PC -XG*EMAGF*A3*PC+ALPHAG*EMAGO* 
@ EMAGF+EMAGF*RX2*C 1 ) / ( RX2*ALPHAG -EMAGO*C 1 -XG*B 3*PC ) 

IF (RRFl .GT. 0.) GO TO 40 
RRF 1=0. 

DPOC5=0 . 

IF (RRF2 .LE. 0.) GO TO 41 
GO TO 42 

40 IF (RRF2 .GT. 0.) GO TO 42 

41 RRF2=0. 

DP0C8=0 . 

42 DPF 1= ( SQRT (PC I**2+WF*PC I*RRF 1)-PCI)/PC 
DPF2=(SQRT(PCI**2+WF*PCI*RRF2) -PCI ) /PC 
GO TO 44 

43 DPFC=(SQRT(PCI**2+WF*PCI*RRF)-PCI)/PC 

44 IF ((DPOC2 .EQ. 0.) .AND. (DPFC .EQ. 0.) .AND. (DPOC5 .EQ. 0.) 
@ .AND. (DPFl .EQ. 0.) .AND. (DPOC8 .EQ. 0.) .AND. (DPF2 

@ .EQ. 0.)) GO TO 27 

C THE LOW AND HIGH FREQUENCY SLOUTIONS ARE SORTED IN 
C ORDER OF INCREASING FUEL PRESSURE DROP RATIO 

IF (NN.EQ. n GO TO 47 
IF (DPFl .EQ. 0.) GO TO 45 
IF (DPF2 .EQ. 0.) GO TO 46 
N=N+1 
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UPL(N)=DPF1 

FL(N)=FREQ 

ACRL(N)=DP0C5 

N=N+1 

UPL(N)=DPF2 
FL(N)=FREQ 
ACRL(N)=DP0C8 
GO TO 27 

45 N=N+1 
UPL(N)=DPF2 
FL(N)=FREQ 
ACRL(N)=DPOC8 
GO TO 27 

46 N=N+1 
UPL(N)=DPF1 
FL(N)=FREQ 
ACRL(N)=DPOC5 
GO TO 27 

47 IF (DPFl .EQ. 0.) GO TO 48 

IF (DPF2 .EQ. 0.) GO TO 49 

N=N+1 

UPH(N)=DPF1 

FH(N)=FREQ 

ACRH(N)=DP0C5 

N=N+1 

UPH(N)=DPF2 
FH(N)=FREQ 
ACRH(N)=DP0C8 
GO TO 27 

48 N=N+1 
UPH(N)=DPF2 
FH(N)=FREQ 
ACRHCN)=DP0C8 
GO TO 27 

49 N=N+1 
UPH(N)=DPF1 
FH(N)=FREQ 
ACRH(N)=DP0C5 
GO TO 27 

50 IF ((FREQ .GT. 300.) .AND. (N .EQ. 0)) JM=1 

N2=N 

IF ((JN .EQ. 1) .AND. (JM .EQ. 1)) GO TO 78 
IF (JM .EQ. 1) GO TO 51 
IF (JN .EQ. 1) GO TO 52 
IM=3 

GO TO 53 

51 IM=1 

GO TO 53 

52 IM=2 

GO TO 56 
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53 


N11=N1-1 
DO 55 IX=1, Nil 
IX1=IX+1 

DO 55 JX=IX1, N1 
IF (UPL(IX)-UPL(JX)) 55,55,54 

54 TEMP=UPL(IX) 

UPL(IX)=UPL(JX) 

UPL(JX)=TEMP 

TEMP=ACRL(IX) 

ACRL(IX)=ACRL(JX) 

ACRL(JX)=TEMP 

TEMP=FL(IX) 

FL(IX)=FL(JX) 

FL(JX)=TEMP 

55 CONTINUE 

IF (IM .EQ. 1) GO TO 59 

56 N22=N2-1 

DO 58 KX=1, N22 
KX1=KX+1 

DO 58 LX=KX1, N2 
IF (UPH(KX)-UPH(LX)) 58, 58, 57 

57 TEMP=UPH(KX) 

UPH(KX)=UPH(LX) 

UPH(LX)=TEMP 

TEMP=ACRH(KX) 

ACRH(KX)=ACRH(LX) 

ACRH(LX)=TEMP 

TEMP=FH(KX) 

FH(KX)=FH(LX) 

FH(LX)=TEMP 

58 CONTINUE 

C THE VARIABLE AREA CASE (JJ=1) CALLS INTERPOLATION TO FIND 
C CRITICAL OX REAL PART AND FREQUENCY AT THE OPERATING FUEL 
C PRESSURE DROP RATIO. THE RESULT IS STORED. 

C THE THROTTLING CASE (JJ=0) CALLS FOR WRITING OUT SOLUTION 

59 IF ((JJ .EQ. 1) .AND. (IM .EQ. 1)) GO TO 60 
IF ((JJ .EQ. 1) .AND. (IM .EQ. 2)) GO TO 62 
IF ((JJ .EQ. 1) .AND. (IM .EQ. 3)) GO TO 64 
GO TO 68 

60 IF ((DPF .LT. UPL(l)) .OR. (DPF .GT. UPL(Nl))) GO TO 79 
LI=LI+1 

NS=N1-1 
DO 61 L=l, NS 

IF (DPF .GT. UPL(L+D) GO TO 61 

ACCL(LI)=ACRL(L)+(DPF-UPL(L))*(ACRL(L+1)-ACRL(L))/(UPL(L+D- 
@ UPL(D) 

FFL(LI)=FL(L)+(DPF-UPL(L))*(FL(L+1) -FL(L) )/ (UPL(L+1) -UPL(L) ) 
UPPL(LI)=DPF 


107 



GO TO 79 

61 CONTINUE 
LI=LI-1 
GO TO 79 

62 IF ((DPF .LT. UPH(l)) .OR. (DPF .GT. UPH(N2))) GO TO 79 
rL=iL+i 

NS=N2-1 
DO 63 L=l, NS 

IF CDPF .GT. UPH(L+D) GO TO 63 

ACCH(IL)=ACRH(L)+(DPF-UPH(L))*(ACRH(L+1)-ACRH(L))/(UPH(L+1) 

@ -UPHCL)) 

FFH(IL)=FH(L)+(DPF-UPH(L))*(FH(L+1)-FH(L))/(UPHCL+1)-UPH(L)) 

UPPH(IL)=DPF 

GO TO 79 

63 CONTINUE 
IL=IL-1 
GO TO 79 

64 IF ((DPF .LT. UPL(l)) .OR. (DPF .GT. UPL(Nl))) GO TO 66 
LI=LI+1 

NS=N1-1 
DO 65 L=l, NS 

IF (DPF .GT. UPL(L+1)) GO TO 65 

ACCL (LI )=ACRL (L)+ (DPF -UPL(L) )* (ACRE (L+1 ) -ACRE (L) )/ (UPL (L+1 ) 

@ -UPL(D) 

FFL(LI)=FL(L)+(DPF-UPL(L))*(FL(L+1)-FL(L))/(UPL(L+1)-UPL(L)) 

UPPL(LI)=DPF 

GO TO 66 

65 CONTINUE 
LI=LI-1 

66 IF ((DPF .LT. UPH(l)) .OR. (DPF .GT. UPH(N2))) GO TO 79 
IL=IL+1 

NS=N2-1 

DO 67 LJ=1, NS 

IF (DPF .GT. UPH(LJ+1)) GO TO 67 

ACCH ( IL)=ACRH (LJ)+ (DPF-UPH (LJ) )* (ACRH (LJ+1) -ACRE (LJ) ) / 

@ (UPH(LJ+1)-UPH(LJ)) 

FFH(IL)=FH(LJ)+(DPF-UPH(LJ))*(FH(LJ+1)-FH(LJ))/(UPH(LJ+1) 

(3 -UPH(LJ)) 

UPPH(IL)=DPF 
GO TO 79 

67 CONTINUE 
IL=IL-1 
GO TO 79 

C WRITE OUT ORDERED SOLUTIONS 

68 WRITE flOUT, 69) 

69 FORMAT (6X, 'FREQUENCY' ,9X, 'DELTAP FUEL/PC' ,6X, 'DELTAP OX/PC' , 
@ 6x , ' FREQUENCY ' , 9x , ' DELTAP FUEL/PC ' ; 9x , ' DELTAP OX/PC ' , / ) 

IF (N1 .GT. 40 .OR. N1 .EQ. 40) ISTEP=40/2 
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IF (N1 .LT. 40) ISTEP=N1 

WRITE (lOUT, 70) (FL(KL), UPL(KL), ACRL(KL), FL(KL+ISTEP) , 
@ UPLCKL+ISTEP) , ACRL(KL+ISTEP), KL=1, ISTEP) 

70 FORMAT (2 (SX,F5 . 1 .2X,E20 . 3 , 2X,E20 . 3) ) 

WRITE (lOUT, 1) 

WRITE (10, 71) (FL(IND), ACRL(IND), UPL(IND), DPO, DPF, 

@ IND=1, 40) 

71 FORMAT (IH ,4X,F5 . 1,2X,E10 .4, 2X,E10 . 4,2X,E10 . 4,2X,E10 . 4) 

C VARIABLE AREA IS SPECIFIED 

78 IF (JJ .EQ.O) GO TO 89 

79 J=2 
K=K+1 
YY=K 

IF (YY .GT. 8.) GO TO 80 
IF (YY .EQ. 1.) GO TO 87 
APF=APF*YY/(YY-1.) 

GO TO 14 

C THE VARIABLE AREA VALUES OF CRITICAL OX REAL PART 
C AND FREQUENCY AT OPERATING FUEL PRESSURE DROP RATIOS ARE 
C CONVERTED TO PROPER VARIABLES FOR ORDERING AND WRITING 

80 IF ((IL .EQ. 0) .AND. (LI .EQ. 0)) GO TO 89 
IF (IL .EQ. 0) GO TO 83 

IF (LI .EQ. 0) GO TO 85 
DO 81 IZ=1, LI 
UPL(IZ)=UPPL(IZ) 

ACRL(IZ)=ACCL(IZ) 

FL(IZ)=FFL(IZ) 

81 CONTINUE 
N1=LI 
IM=3 
JJ=0 

DO 82 IZ=1, IL 
UPH(IZ)=UPPH(IZ) 

ACRH(IZ)=ACCH(IZ) 

FH(IZ)=FFH(IZ) 

82 CONTINUE 
N2=IL 

GO TO 53 

83 DO 84 IZ=1, LI 

UPL(IZ)=UPPL(IZ) 

ACRL(IZ)=ACCL(IZ) 

FL(IZ)=FFL(IZ) 

84 CONTINUE 
N1=LI 
IM=1 
JJ=d 
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85 


GO TO 53 
DO 86 IZ=1, IL 
UPH(IZ)=UPPH(IZ) 

ACRH(IZ)=ACCH(IZ) 

FH(IZ)=FFH(IZ) 

86 CONTINUE 
N2=IL 
IM=2 
JJ=0 

GO TO 56 

87 APF=.25*APF 
GO TO 14 

88 J=1 

89 10=10+1 

90 CONTINUE 
WRITE (5, 92) 

92 FORMAT (IH , 'ANOTHER STUDY?' ,//,3x, ' 1 YES',/,3x,'2 NO') 

READ (5, *) INPUT 
IF (INPUT -EQ. 1) THEN 
WRITE (lOUT, 94) 

94 FORMAT (IHl) 

GO TO 15 
END IF 
STOP 
END 


FUNCTION TM(COEFF) 

C THIS MAPS THE EMPIRICAL MIXING CURVE 
DIMENSION S(7), V(7) 

DATA S(l)/.5/,S(2)/l./,S(3)/1.5/,S(4)/1.8/,S(5)/2./,S(6) 

(§ 11.11 

DATA V(l)/.15/,V(2)/.35/,V(3)/.7/,V(4)/1.075/,V(5)/1.4/, 
@ V(6)/1.9/ 

DO 1 1=1, 5 

IF ((S(I) .LE. COEFF) .AND. (COEFF .LT. S(I+1))) GO TO 2 

1 CONTINUE 

2 TM=1.E-3*(V(I)+(C0EFF-S(I))*(V(I+1)-V(I))/(S(I+1)-S(I))) 
RETURN 

END 
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SUBROUTINE GAM (TC, GAMMA) 

C THIS SUBROUTINE COMPUTES THE AVERAGE SPECIFIC HEAT RATIO OF THE 
C PRODUCTS USING THE EQUATIONS FOUND IN VAN WYLEN. 

THETA=TC/100. 

R=1.987 

CPH2=13. 505-167. 9 6*THETA** ( - 0 . 7 5 ) +2 7 8 . 44*THETA** ( - 1 ) - 1 34 . 0 1 
@ *THETA**(-1.5) 

CPH20=34 . 19 -43 . 868*THETA^>’'KO . 25 )+19 . 778*THETA**(0 . 5 ) -0 . 88407 
<a *THETA 

CPAVG=(5./20.)*CPH2 + (1 . /20 . )*CPH20 
GAMMA=1 . / ( 1 . -R/CPAVG) 

RETURN 

END 
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